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I.  INTRODUCTION 


The  theory  of  integral  equations  is  a  rich  and  complex  theory  that  touches 
every  branch  of  science  and  engineering.  Despite  its  prevalence,  few  attempt  to 
master  the  theory  of  integral  equations.  The  history  of  integral  equations  dates 
back  to  the  early  nineteenth  century  when  the  profound  mathematical  insights  of 
Newton  and  Leibniz  were  being  used  in  conjunction  with  rapid  technological 
advances  to  describe  physical  phenomena  that  had  defied  scholars  for  millennia. 
While  most  texts  offer  a  brief  history  of  integral  equations  in  the  introduction, 
Lonseth  [1]  provides  a  paper  that  is  rich  with  history  and  mathematical  rigor. 

A  fundamental  problem  in  mechanical  engineering  is  the  efficient 
prediction  for  transient  responses  of  components  and  structures,  especially  after 
modifications  have  been  made.  The  design/analysis  cycle  often  requires 
repeated  analyses,  as  the  designer  looks  for  an  optimal  design.  Structural 
synthesis  is  the  analysis  of  the  dynamic  response  of  a  system  when  either 
subsystems  are  combined  (substructure  coupling),  or  modifications  are  made  to 
subsystems  (structural  modification)  [2],  Structural  modification  can  include 
changes  in  properties  such  as  mass,  stiffness  or  damping.  Substructure  coupling 
can  occur  for  many  reasons.  It  can  be  a  natural  product  of  the  system  itself  as  a 
result  of  factors  such  as  material  differences  or  complex  geometries. 
Substructure  coupling  is  frequently  used  when  multiple  internal  or  external 
organizations  individually  contribute  subsystems  to  a  governing  system  design. 
The  goal  when  conducting  this  analysis  is  to  maximize  computational  efficiency, 
a  key  requirement  to  remain  competitive  on  a  global  scale. 

There  are  many  methods  for  structural  synthesis  and  the  method  this 
paper  focuses  on  is  the  integral  equation  formulation  for  transient  structural 
synthesis  as  formulated  by  Gordis  in  [2],  This  formulation  is  beneficial  because  it 
eliminates  the  process  of  reconstructing  and  resolving  large  systems  of 
equations  after  a  modification  to  the  structure  has  been  made.  Also,  the 

response  to  the  modified  system  is  based  on  a  pre-modified,  or  “baseline” 
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system  response  usually  calculated  from  a  pre-existing  model.  This  bypasses  the 
process  of  reassembling  and  solving  for  the  solution  to  the  modified  system. 
Furthermore,  synthesis  can  reduce  the  scope  of  the  analysis  from  the  total 
structure,  to  the  portions  of  the  structure  in  which  the  modification  has  been 
made.  Responses  can  be  calculated  from  any  portion  of  the  structure  that  has 
not  been  modified  but  is  of  interest. 

The  integral  equation  formulation  to  structural  synthesis  is  based  on  the 
convolution  integral  and  results  in  a  Volterra  integral  equation  (VIE)  of  the  second 
kind.  Integral  equations  usually  cannot  be  solved  analytically  and  one  must  resort 
to  numerical  schemes.  Since  integral  equations  can  be  re-written  in  the  form  of 
differential  equations,  it  implies  that  same  numerical  methods  used  to  solve 
differential  equations  can  be  used  to  solve  integral  equations.  Accuracy  varies 
with  the  various  methods,  but  in  general,  the  higher  the  order  of  accuracy,  the 
more  complex  the  method.  The  approach  this  paper  uses  to  solve  the  integral 
formulation  for  transient  structural  synthesis  is  an  adaptive  method  based  on  the 
trapezoidal  rule  formulated  by  Gordis  and  Neta  [3], 

This  method  computes  the  whole  integral  as  the  sum  of  integrals 
evaluated  on  sub-intervals  between  the  limits  of  integration.  The  trapezoidal  rule 
is  applied  to  each  of  these  sub-integrals  and  when  summed  together  forms  the 
coarse  solution.  A  fine  solution  is  created  by  dividing  the  final  sub-interval  into 
two  equal  intervals,  raising  the  accuracy  level  of  the  final  interval  and  therefore 
the  overall  solution.  This  is  necessary  only  for  the  final  interval  because  all  other 
previously  calculated  solutions  are  assumed  to  have  been  calculated  to  within  an 
acceptable  level  of  accuracy.  The  final  interval  is  the  only  interval  that  has  the 
value  from  the  coarse  solution  in  which  the  accuracy  is  in  question.  If  the  error 
between  the  fine  and  course  solution  meets  the  specified  tolerance,  the  fine 
solution  is  stored  as  the  coarse  solution,  the  time-step  is  doubled,  and  the 
process  repeated  for  the  next  interval.  If  the  error  is  greater  than  the  tolerance, 
the  solution  at  the  current  point  is  deleted,  the  time-step  is  halved,  and  the 
process  repeated  from  the  last  point  of  acceptable  accuracy. 
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The  adaptive  trapezoidal  method  is  beneficial  because  the  level  of 
accuracy  can  be  established  and  the  time-step  then  adjusted  as  necessary  to 
meet  the  desired  level  of  accuracy.  The  only  unknown  that  must  be  determined 
by  this  method  is  the  intermediate  function  value  in  the  final  interval  which  is  a 
straightforward  calculation.  Higher  order  methods  require  a  larger  number  of 
intermediate  values  in  the  final  interval  and  therefore  a  higher  number  of 
intermediate  function  values  that  need  to  be  solved.  The  higher  level  of  accuracy 
normally  obtained  by  higher  order  methods  is  obtained  by  the  adaptive 
trapezoidal  rule  when  the  time-step  is  sufficiently  small,  or  when  the  time-step  is 
reduced  in  order  to  meet  the  set  error  tolerance.  In  [4],  Linz  gives  a  detailed 
analysis  of  the  trapezoidal  rule  and  other  higher  order  methods,  along  with  the 
truncation  error  of  each  method. 

Initially  the  goal  was  to  verify  that  the  adaptive  trapezoidal  formulation 
would  be  a  robust  alternative  to  solving  synthesis  problems.  After  successfully 
duplicating  the  results  of  [3],  the  adaptive  method  was  successfully  applied  to  the 
synthesis  of  a  single  degree-of-freedom  (SDOF)  mass-spring  system  with  a 
stiffness  modification,  and  a  two  DOF  mass-spring  system  with  a  stiffness 
modification.  Both  of  these  problems  were  subjected  to  step  and  periodic  forcing 
functions.  However,  when  applying  this  method  to  an  aluminum  cantilever  beam 
subjected  to  a  step  forcing  function  with  an  elemental  stiffness  change,  the 
method  exhibited  conditional  stability.  Upon  further  investigation,  this  conditional 
stability  is  attributed  to  the  magnitude  of  the  stiffness  modification. 

Stability  analysis  of  integral  equations  is  complex  and  is  not  the  focus  of 
this  paper.  For  those  interested,  the  issue  is  addressed  in  great  detail  in  [5],  [6], 
and  [7],  In  general,  integral  equations  with  large  kernels  are  not  unconditionally 
stable  [7],  In  structures  with  modifications  of  sufficient  magnitude,  the  problem 
becomes  conditionally  stable.  Since  the  modification  is  a  component  of  the 
kernel;  therefore  the  magnitude  of  the  kernel  can  become  large  enough  to  drive 
the  problem  from  unconditional  stability  to  conditional  stability. 
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In  order  to  overcome  conditional  stability,  the  derivative  of  the  synthesis 
equation  was  taken,  producing  a  Volterra  integro-differential  equation  (VIDE)  of 
the  second  kind.  Stability  of  the  original  synthesis  equation  and  the  motive  for 
taking  the  derivative  will  be  discussed  in  greater  detail  in  the  Theory  section  of 
the  paper.  Though  taking  the  derivative  complicates  the  adaptive  algorithms, 
taking  the  derivative  of  the  synthesis  equations  has  three  benefits: 

1 )  Provides  unconditional  stability  to  the  synthesis  equation, 

2)  Provides  velocity  information  that  was  previously  unavailable,  and 

3)  The  solution  to  the  original  synthesis  equation  is  inherent  in  the  algorithm 
when  solving  for  the  derivative. 
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II.  THEORY 


The  main  focus  of  this  analysis  will  be  single  degree-of-freedom  and 
multiple  degree-of-freedom  (MDOF)  systems  that  have  undergone  structural 
modifications;  specifically  stiffness  modifications.  Also,  this  analysis  is  strictly 
limited  to  linear  models.  The  theory  of  transient  structural  modification,  and  the 
associated  governing  VIE,  is  based  on  the  convolution  integral,  or  Duhamel’s 
integral.  The  convolution  integral  is  commonly  used  when  analyzing  a  system 
response  to  an  external  excitation.  For  structural  modifications  that  are  of 
sufficient  magnitude  this  VIE  is  conditionally  stable  relative  to  the  magnitude  of 
the  stiffness  modification,  but  it  will  be  shown  that  converting  the  VIE  to  a 
Volterra  integro-differential  equation  provides  unconditional  stability;  allowing  the 
use  of  the  adaptive  trapezoidal  rule  to  solve  the  system. 

A.  INTEGRAL  EQUATIONS 

Prior  to  discussing  the  integral  formulation  of  structural  synthesis,  it  is 
essential  to  understand  the  basic  structure  of  integral  equations  in  general  terms. 
The  following  discussion  can  be  found  in  greater  detail  in  any  text  on  integral 
equations,  but  this  discussion  follows  closely  from  [8], 

An  integral  equation  is  expressed  generally  as: 

b(x) 

h(x)u(x)  =  f(x)+  J  K(x,4)u(4)dg  (1) 

0 

When  b(x)  is  constant,  the  equation  is  classified  as  a  Fredholm  equation  and  is 
not  the  focus  of  this  paper.  When  b(x)  =  x ,  the  equation  is  classified  as  a  Volterra 
equation.  Furthermore,  if  b(x)  =  0,  then  both  Fredholm  and  Volterra  equations 

are  of  the  first  kind.  If  h(x)  =  c,  where  c  is  a  constant,  then  both  Fredholm  and 

Volterra  equations  are  of  the  second  kind.  As  previously  mentioned,  the 
convolution  integral  used  to  formulate  structural  synthesis  takes  the  form: 
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(2) 


X 

u(x)  =  f(x)  +  Jtf  (x,  g)u(g)  dg 
0 

a  Volterra  equation  of  the  second  kind  where  f(x)  is  a  known  function,  K{x,£) 
is  the  known  kernel  function,  and  u  is  the  function  to  be  solved  for. 

B.  THE  CONVOLUTION  INTEGRAL  (DUHAMEL’S) 

The  convolution  integral,  or  Duhamel’s  integral  expression  can  be  used  to 
describe  the  response  of  a  system  to  an  external  excitation  and  is  based  on  the 
principle  of  superposition;  which  is  only  applicable  to  linear  systems  [9],  The 
Duhamel’s  integral  expression  in  terms  of  dynamic  response  is  as  follows: 

MO} = MM, +\[h(t  -  r)]  M  (T)}dr  (3) 

o 

where  jx(r)}  is  a  nx  1  vector  of  the  total  system  response  (where  n  is  the 
number  of  DOF  of  the  system),  M  (0}  's  a  nx 1  vect°r  of  the  free  response  of 
the  system,  [h(t- r)J  is  an  nxn  impulse  response  function  (IRF)  matrix,  and 

|/c(r)|  is  a  nx  1  vector  of  external  forces  (external  loading)  applied  to  the 
system. 

In  terms  of  \xh  (r)] ,  the  solution  of  the  free  response  will  vary  depending  on 

whether  the  system  is  underdamped,  critically  damped,  or  overdamped.  The 
equation-of-motion  (EOM)  for  a  linearly  viscous  damped  SDOF  system  is: 

x  +  2i;(onx  +  afnx  =  0  (4) 

and  assuming  the  solution  to  the  ordinary  differential  equation  is  an  exponential, 
the  characteristic  polynomial  can  be  solved  and  the  eigenvalues  determined  from 
the  following: 
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(5) 


It  can  be  seen  that  the  magnitude  of  the  damping  ratio  (<^)  determines  the 

character  of  the  free  response  EOM.  There  are  three  possible  cases  for  the 
damping  ratio:  1 )  0<  ^  <1  (underdamped),  2)  C,  =1  (critically  damped),  and  3)  C, 
>1  (overdamped).  Substituting  the  eigenvalues  into  the  exponential  solution  and 
applying  initial  conditions  gives  the  solution  for  the  three  cases  of  damping. 

x(t)  =  e  x0  cos ( wdt )  +  — - —  sin ( wdt )  ( underdamped)  (6) 

\  ) 

x(t)  =  [x0  +  (x0  +  (Onx0 )/]  c  <’<°n>  ( critically  damped)  (7) 

x (t )  =  c  x0cosh ( co  t  j  +  io  +  sinh (cojt^  (o ver damped)  ( 8) 

In  any  case,  following  the  assumptions  thatx0  =  x0  =  0  ,  this  leads  to  {x/7  (r)}  ={0} 
and  therefore  the  convolution  integral  becomes: 

t 

x(t)  =  \h(t-T)f(r)dT  (9) 

0 

In  general,  structural  systems  are  underdamped  and  therefore  will  be  the  focus  of 
this  paper. 

C.  IMPULSE  RESPONSE  FUNCTION 

It  is  important  to  understand  the  impulse  response  function  matrix 
because  it  is  the  kernel  of  the  VIE  formulation  of  the  synthesis  equation;  and  is 
the  component  of  the  VIE  that  determines  stability  of  the  system. 

The  EOM  of  a  MDOF,  proportionally  damped  system  with  an  external 
excitation  applied  to  a  single  DOF  is  as  follows: 

[M]{x}  +  [c]|x}  +  [^]|x}  =  {p(r)}  (10) 

Often  this  system  of  equations  are  coupled,  in  other  words,  there  are  not  enough 
independent  equations  to  solve  for  the  entire  system.  A  transformation  is 
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introduced  that  takes  the  EOM  from  physical  coordinates  to  modal  coordinates. 
This  transformation  uncouples  the  equations  allowing  for  the  appropriate  number 
of  independent  equations  to  solve  for  the  system.  The  transformation  that  takes 
the  EOM  from  physical  coordinates  to  modal  coordinates  is  as  follows: 

W  =  WW  (n) 

Substituting  Equation  (11)  into  Equation  (10),  and  multiplying  the  EOM  by[®]7, 
gives  the  following  EOM  in  modal  coordinates: 

[Afl{«}+[c]{4}+M{«}=[*f{i>(0}  m 

where,  [A1]  =  [0]r[M][®],[C]  =  [®]r[C][®],  and [/C]  =  [<f>]7  [^][<I>].  Assuming  that 

the  system  is  mass-normalized,  and  the  external  excitation  acts  at  a  single  DOF, 
Equation  (12)  becomes: 

O' 

0 

{<?}+  2^yn  {q}+  co]  {q}  =  [of «  1  >  p(t)  (13) 

0 

0 

where  1  is  at  the  jth  DOF.  Equation  (13)  then  reduces  to: 

{<?}  +  2  <Za)n  {4}+  co]  (?}  =  {^}p(0  (W 

Since  the  system  in  modal  coordinates  is  uncoupled,  each  equation  can  be 
solved  independently.  The  jth  modal  EOM  can  be  written  as: 

q,  +  2  CiO)^  +  (o]qt  =  (j)fp  (t)  (15) 

In  order  to  determine  the  /"'impulse  response  function  Equation  (15) 
needs  to  be  re-written  as: 
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(a) 

(b) 


(16) 


\p(t)</>o\0<t<td 

jo  ,td<t 


where  td  is  the  time  the  impulse  is  applied.  Recall  the  impulse  is  defined  as: 

ld 

I  =  J  p(t)dt  (17) 

o 


Taking  the  integral  of  Equation  16(a)  from  0  to  td : 

rd  rd  ld  ld 

J  +  J  2Cicon  q(t)dt  +  J  orn  q(t)dt  =  J (jy'd p{t)dt 


(18) 


Given  the  initial  rest  conditions  of  q(0)  =  0  and  4(0)  =  0,  Equation  (18)  becomes: 


[q{td  )  -  4(0)]  +  2  4>,,;  \q(td  )  -  q(0)]  +  (o\  qavgtd  =  I(^ 


(19) 


Applying  the  initial  conditions  and  understanding  that  qUVK  is  very  small  while  the 
impulse  is  applied,  the  last  term  can  be  ignored  and  Equation  (19)  becomes: 

)(td )  +  2  q(td )  =  Hf  (20) 

In  the  short  duration  the  impulse  is  applied,  g(0)  =  0and  4(0)  =  which  become 

the  initial  conditions  when  solving  for  Equation  16(b).  Solving  Equation  16(b) 
yields: 

q(t)=  c  sinico p)  (21) 

™di  V  '  7 

Letting  the  impulse,  1  =  1,  the  modal  IRF  is: 

(h(i) 

H(t)  =  ^—ei‘a>',‘lsin(colt )  (22) 

w  oj,  v  '  ’ 


To  transform  the  modal  IRF  to  the  physical  IRF,  Equation  (11)  is 
substituted  into  Equation  (21).  Applying  this  transformation  to  the  system  of 
modal  IRF  gives  the  physical  IRF  represented  in  matrix  form  as: 


[MOM®] 


(23) 
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(24) 


From  Equation  (23)  it  follows  that  the  IRF  matrix  entries  take  the  form  of: 

hijit  -t)  =  —e^C0"(t-T)  sin (o)d(t  -  r)) 

p= i  md 

It  is  important  to  keep  in  mind  that  the  physical  parameters  of  the  IRF 
matrix  {C,,(on,  andry^)  are  determined  by  the  baseline  model.  It  should  also  be 
pointed  out  that  Equation  (24)  applies  to  an  underdamped  MDOF  system  (<^<1) 
with  no  rigid  body  modes;  a  reasonable  assumption  for  a  cantilevered  beam. 

D.  INTEGRAL  EQUATION  FORMULATION  FOR  STRUCTURAL 
MODIFICATION 

When  a  structural  modification  has  been  made,  the  convolution  integral 
will  model  the  modified  response  jx*(r)j: 

{x*(t)}=\[h(t-r)]{f(T)}dT  (25) 

0 

where  {/(r)}  is  the  sum  of  the  externally  applied  excitations  and  the  new  force 
imposed  on  the  system  as  a  result  of  the  system  modification  and  is  as  follows: 

{/«}={/'M}-{/‘«}  (26) 

where  j/e(r)}  's  the  vector  of  externally  applied  excitations,  and  {/*(?)}  is  the 
change  in  force  due  to  structural  modifications  and  can  be  written  as: 

{/(r)}  =  {[AM]{x(r)}  +  [AC]{x(r)}  +  [A^]{x(r)^  (27) 

Substituting  Equation  (21)  and  Equation  (22)  into  Equation  (20),  and 

assuming  there  are  no  mass  or  damping  modifications,  the  modified  convolution 
integral  becomes: 

K  (0} = J  [ h  -  o]  [{/'  ( r)}  ~  [■ K  o r)}] dr  (28) 

0 

{jc*(0}  =  }  [ h  (t  -  T)]{fe  (r)}  dz  -  J  [ h  ( t  -  r)]  [A^]  {x*  (r)}  dr  (29) 

0  0 
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Recognizing  that  the  first  term  is  Equation  (9),  the  baseline  response  of  the 
unmodified  system,  Equation  (23)  can  be  reduced  to: 

|V(r)}  =  |x(r)}-J[/7(r-r)][A^]|x*(r)|dr  (30) 

0 

This  equation  takes  the  form  of  Equation  (2),  a  VIE  of  the  second  kind. 
Recall  that  the  unmodified  system  parameters  are  used  to  construct  the  IRF 
matrix  and  therefore  it  is  clear  that  the  modified  system  response  only  requires 
the  baseline  system  response  and  its  parameters,  and  the  structural  modification. 
Though  this  paper  focuses  on  strictly  stiffness  modifications,  the  analysis  can  be 
extended  to  mass  and  damping  modifications,  or  any  combination  of  structural 
modifications.  This  extends  the  analysis  to  solving  integral  equations  with  first 
and  second  order  terms. 


E.  STABILITY  OF  THE  INTEGRAL  EQUATION  FORMULATION  FOR 

STRUCTURAL  SYNTHESIS 

In  order  to  understand  the  stability  of  the  integral  equation  formulation  for 
structural  synthesis,  one  only  needs  to  study  the  coarse  solution  derivation  for 
the  adaptive  method  when  applied  to  the  structural  synthesis  Equation  (25). 
Since  this  is  a  VIE  of  the  second  kind,  the  process  outlined  by  Gordis  and  Neta 
[3]  can  be  followed.  Since  this  derivation  is  not  trivial,  it  has  been  included  as 
Appendix  A. 

Starting  with  Equation  (25): 

K  (0}  -  MO}  - ]  [h(f  ~ 0] MM  (0} dr 

o 


and  following  the  steps  as  outlined  by  [3],  the  coarse  solution  of  Equation  (25)  is: 

[l]  +  l si  Wjj  [ K } ;  =  "Z T +  Si ) [h]..  [A/l]{x* } .  +  {x}  .  (31) 
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where  8  is  the  time-step  for  a  given  interval.  Also  from  Equation  (19),  it  is  seen 
that  [h]  =0  (when  sin(r-r)  =  0  ).  Equation  (26)  now  becomes: 
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(32) 


{xlj=-Tl(S^+S‘)[h]ji[AK]{xli+{x}j 

7=0  Z 

Equations  (26)  and  (27)  give  the  basis  for  the  discussion  of  stability.  The 
structure  used  in  this  analysis  is  an  aluminum  cantilevered  beam.  The  beam  is 
developed  as  a  finite  element  model  in  which  a  stiffness  modification  can  be 
made  to  any  of  the  elements.  An  excellent  discussion  concerning  the 
construction  of  a  finite  element  model  for  a  cantilevered  beam  is  given  by  Kwon 
and  Bang  [10].  The  key  detail  relevant  to  the  discussion  of  stability  is  the  element 
stiffness  matrix.  As  shown  in  [10],  the  element  stiffness  matrix  is  as  follows: 


"  12 

61 

-12 

61 

El 

61 

4 12 

-61 

2 12 

Z3 

-12 

-61 

12 

-61 

61 

2l2 

-61 

4 12 

The  element  stiffness  modification  to  the  structure  can  be  represented  as: 

[AK]=7,[r]  (34) 

where  q  is  the  percent  change  of  the  element  stiffness.  The  parameters  of  the 
aluminum  beam  are:  Young’s  Modulus (£)  =  69  GPa  (lO 6  psi),  length (/)  =  3.05  m 

(120  in)  ,  depth (d)  =  0.3  m  (12  in)  ,  and  width  (w)  =  0.1  m  (4  in).  This  results  in  a 

moment-of-inertia  (/)  of. 0.01m3  (576  m3).  When  substituting  the  values  into 

Equation  (28),  the  magnitude  of  the  stiffness  change  is  on  the  order  of  M08 . 

Understanding  that  the  coarse  solution  for  the  integral  formulation  for 
structural  synthesis  is  governed  by  Equation  (27),  the  kernel  is  then  defined  as: 

a*) 

The  stiffness  modification  is  a  fixed  value  and  can  be  of  sufficient  magnitude  to 
make  the  VIE  conditionally  stable  for  the  given  time-step.  The  only  parameter 

that  can  be  adjusted  is  choice  of  the  time-step;  affecting  both  (Si  +  4 1)  and  [h] 
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terms  sufficiently  to  return  the  VIE  to  stability.  This  was  observed  running  the 
code  that  modelled  the  cantilevered  beam.  If  the  time-step  was  not  sufficiently 
small  for  a  given  stiffness  modification,  then  the  solution  was  unstable.  A 
sufficiently  small  time-step  would  give  a  stable  and  accurate  solution  but  will 
increase  the  computation  time. 

Trying  to  determine  a  sufficiently  small  time-step  for  stability  is  not 
effective,  especially  when  considering  an  adaptive  time-step  method.  Returning 
to  Equation  (26),  it  is  important  to  realize  that  the  left  hand  IRF  function  matrix 
will  always  be  zero.  A  non-zero  left  hand  IRF  function  would  effectively  normalize 
Equation  (26)  and  maintain  unconditional  stability  of  the  equation.  In  order  to 
obtain  a  non-zero  term  on  the  left-hand  side,  the  logical  step  is  to  take  the 
derivative  of  the  synthesis  equation. 

F.  DERIVATIVE  OF  THE  SYNTHESIS  EQUATION 

The  goal  of  taking  the  derivative  of  the  synthesis  equation  is  to  ensure  that 
the  kernel  does  not  equal  zero  when  (r  =  r).  The  full  derivation  of  the  derivative 

of  |i*|  is  shown  in  Appendix  B,  but  the  resulting  equation  is: 

r 

=  [A/f]{x’"(r)j dr  (36) 

0 

a  VIDE.  Now  the  entries  of  the  derivative  of  the  IRF  matrix  take  the  form: 

K,  « - *■>  =  S  «-!■))+ * cos(n  (f -  r)))  (37) 

P= 1 

When  r  =  t  ,  Equation  (32)  reduces  to: 

4<o) =2>,v;  (38) 

p= i 

Or  in  terms  of  the  system: 

[/i(0)]  =  [O][of 
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(39) 


Applying  the  process  outlined  in  [3],  the  coarse  solution  of  Equation  (31) 


is: 


+  «[*]  M{*}, -\s,[h\  [4A:](f  +  {i-}M))  (40) 


--S] 

4  7 


h]\AK]{x 


1  - 

X  }  -~Sj 
j- 1  2  7 


2'LI^K 


x  )  +\x 
0 


The  full  derivation  of  Equation  (35)  can  be  found  in  Appendix  C.  There  are  a  few 
points  to  notice: 

1)  Equation  (36)  is  substantially  more  complicated  than  Equation  (26). 

2)  Equation  (36)  contains  both  jx*|  and  |x*|  terms. 


3)  The  left-hand  side  of  Equation  (36)  contains  a  non-zero  derivative  of  the 
IRF  matrix.  Substituting  Equation  (34)  into  the  left-hand  side  of  Equation 
(35),  the  left-hand  side  of  Equation  (35)  can  be  re-written  as: 


U  +  -6] 

4  1 


Ji]  [ajcDK), -{*•}, m 


The  first  challenge  is  to  determine  if  Equation  (35)  can  be  solved  in  its 
current  form  with  both  |x*j  and  jx*j  terms  present.  Using  the  trapezoidal  rule, 

|x*|  can  be  written  as  follows: 


(42) 


Due  to  the  iterative  nature  of  the  algorithm,  the  current  modified  displacement, 
jx*j  ,  can  be  solved  as  soon  as  the  velocity,  {**}  .,  for  the  given  time  step  is 

solved.  Furthermore,  the  previously  calculated  modified  displacement,  |x*|  , 

and  previously  calculated  modified  velocity,  jx*}  ^ ,  can  be  used  to  calculate  the 

current  modified  displacement.  This  becomes  clear  after  a  few  iterations: 
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FV!*T-.+f  ({iT-.+(*"},) 


This  highlights  one  of  the  benefits  of  evaluating  the  derivative  of  the 
synthesis  equation.  Even  though  the  resulting  algorithm  is  far  more  complicated 
than  the  original  synthesis  equation,  not  only  is  the  velocity  obtained,  the 
modified  displacement  is  also  obtained. 

The  next  question  to  be  answered  is  to  whether  the  non-zero  derivative  of 
the  IRF  matrix  on  the  left-hand  side  returns  the  VIDE  to  unconditional  stability. 

The  term  [[/]+ ^;  [*]yj  [A K]j  is  analogous  to  a  scalar  in  a  SDOF  problem. 

Performing  the  appropriate  operations  to  solve  for  ji*j  is  also  analogous  to 

dividing  the  equation  by  the  scalar  in  a  SDOF  problem.  This  normalizes  the 
equation  at  each  time-step,  essentially  reducing  the  effect  that  the  stiffness 
modification  has  on  the  system  and  maintaining  unconditional  stability. 

Determining  the  stability  of  integral  equations  is  not  a  trivial  task  and  is 
usually  addressed  on  a  case-by-case  basis.  In  order  to  determine  whether  the 
derivative  of  the  synthesis  equation  was  unconditionally  stable,  a  range  of  time- 
steps  was  tested.  If  the  system  was  still  conditionally  stable,  at  a  large  time-step 
the  system  would  produce  an  unstable  result.  Running  the  program  with  an  initial 
time-step  of  1.0x10  3 sec  produced  a  stable,  but  inaccurate  solution  that  closely 
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modelled  the  unmodified  system  response.  The  step-size  was  reduced  by  a 
magnitude  of  10  until  the  time-step  of  MO  5 sec  was  reached.  At  this  time-step, 
the  solution  closely  followed  the  exact  solution  to  the  modified  system  response. 
Not  only  is  the  system  unconditionally  stable,  but  it  produces  accurate  results  for 
the  modified  system  with  the  appropriate  time-step. 
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III.  MODELS  AND  RESULTS 


A.  GENERAL  PROCESS 


The  following  models  have  been  solved  using  the  adaptive  solution 
procedure: 


1)  A  SDOF  mass-spring  system  with  an  externally  applied  step  excitation, 
subjected  to  a  stiffness  modification. 


mass  ( m )  =  1  kg  (2.2  lb),  spring  constant  (k)  =  100 N  /  in 


|^6.9 


lb' 

fi) 


2)  A  SDOF  mass-spring  system  with  an  externally  applied  periodic  excitation, 
subjected  to  a  stiffness  modification. 


mass  (m)  =  1  kg  (2.2  lb),  spring  constant  (k)  =  10(W/  m 


l  ft) 


3)  A  two  DOF  mass  spring  system  with  an  externally  applied  step  excitation, 
subjected  to  a  stiffness  modification  to  the  second  spring. 


mass  (m)  =  1  kg  (2.2  lb),  spring  constant  (k)  =  1002//  m 


l  ft) 


4)  A  generalized  MDOF  cantilevered  aluminum  beam  with  an  externally 
applied  step  excitation,  subjected  to  a  stiffness  modification  to  an  arbitrary 
beam  element.  Recall  this  is  the  model  in  which  conditional  stability  has 
been  observed. 


5)  The  same  model  as  4,  but  the  derivative  of  the  synthesis  equation  has 
been  analyzed  in  order  to  determine  if  unconditional  stability  is  observed. 

The  program  to  the  adaptive  method  has  three  distinct  phases: 

1)  Programming  algorithms  for  the  coarse,  half,  and  fine  approximations  as 
derived  from  Gordis  and  Neta  [3], 

2)  Determining  the  error  between  the  coarse  and  fine  approximation. 

3)  Halving  or  doubling  the  time-step  as  determined  by  the  error. 
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Since  the  algorithms  of  [3]  were  derived  for  the  generalized  VIE  of  the 
second  kind,  the  algorithms  had  to  be  altered  for  the  specific  model.  In  order  to 
ensure  that  the  algorithms  were  correct,  a  non-adaptive  approach  was 
implemented.  Since  the  exact  solution  for  all  models  is  known,  the  results  of  the 
course,  half,  and  fine  approximations  were  calculated  and  individually  compared 
against  the  exact  solution.  Once  these  algorithms  were  properly  programmed, 
then  the  error  assessment  and  adaptive  time-step  adjustment  phases  were 
implemented.  These  algorithms  did  not  significantly  deviate  from  model-to-model. 

As  previously  mentioned,  the  adaptive  solution  for  the  cantilever  beam 
exhibited  conditional  stability  while  evaluating  the  coarse  approximation.  At  this 
point  in  the  analysis,  the  coarse  approximation  was  not  yet  adaptive  and  is 
simply  a  standard  trapezoid  rule.  Once  conditional  stability  had  been  exhibited  in 
the  coarse  approximation,  the  half  and  fine  approximations  were  not  developed 
and  the  derivative  of  the  synthesis  equation  was  pursued.  In  developing  the 
coarse  approximations  to  the  derivative,  it  became  clear  that  though  this  had 
produced  an  unconditionally  stable  and  accurate  solution,  the  algorithm  is  limited 
by  computational  efficiency.  As  a  result,  the  half  and  fine  approximations  were 
not  developed  for  this  model  as  well. 

B.  SDOF  MASS-SPRING  SYSTEM  WITH  EXTERNALLY  APPLIED  STEP 

FUNCTION  SUBJECTED  TO  A  STIFFNESS  MODIFICATION 

The  physical  representation  of  this  model  is  given  in  Figure  1.  The  stiffness 


Spring 

Modification 


Figure  1 .  SDOF  mass-spring  with  stiffness  modification 
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modification  is  varied  in  the  following  intervals:  0%,  10%,  25%,  50%,  75%,  and 
100%  at  an  initial  time-step  of  0.01  second,  and  a  maximum  error  between 
coarse  and  fine  approximations  of  M0-4.  Table  1  gives  the  accuracy,  total 
number  of  calculations,  and  time  to  complete  for  each  stiffness  modification.  The 
plots  of  modified  displacement  and  step  size  for  10%,  50%  and  100%  have  been 
included  to  show  the  variance  in  system  response  as  the  stiffness  modification  is 
increased. 


Stiffness 

Greatest 

Least 

Total  number 

Time  to 

Modification  (%) 

Error  (%) 

Error  (%) 

of  calculations 

complete(sec) 

0 

0.9 

0.0 

7 

0.16 

10 

1.3 

0.0 

902 

3.7 

25 

2.6 

0.0 

1184 

5.0 

50 

4.5 

0.0 

1657 

7.1 

75 

6.2 

0.0 

1803 

7.9 

100 

7.6 

0.0 

1913 

8.5 

Table  1 .  Response  to  stiffness  modifications  for  SDOF  mass-spring 
system  subjected  to  externally  applied  step  excitation. 


x*(t)  vs.  Time 


Figure  2.  Modified  response  for  a  SDOF  mass-spring  with  externally 
applied  step  excitation  subjected  to  10%  stiffness  modification. 
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Step  Size  vs.  Step  Number 


Figure  3.  Step  size  for  a  SDOF  mass-spring  with  externally  applied  step 
excitation  subjected  to  10%  stiffness  modification. 


x*(t)  vs.  Time 


Figure  4.  Modified  response  for  a  SDOF  mass-spring  with  externally 
applied  step  excitation  subjected  to  50%  stiffness  modification. 
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x*(t)  (modified  displacement) 


Step  Size  vs.  Step  Number 


Figure  5.  Step  size  for  a  SDOF  mass-spring  with  externally  applied 
sinusoidal  excitation  subjected  to  50%  stiffness  modification. 


Figure  6.  Modified  response  for  a  SDOF  mass-spring  with  externally 
applied  step  excitation  subjected  to  100%  stiffness  modification. 
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Step  Size  vs.  Step  Number 


Figure  7.  Step  size  for  a  SDOF  mass-spring  with  externally  applied 
sinusoidal  excitation  subjected  to  50%  stiffness  modification. 

Comparing  Figures  2,  4,  and  6  show  that  the  lower  the  stiffness 
modification,  the  more  accurate  the  approximation  to  the  exact  solution.  This  is 
also  evident  in  the  maximum  error  given  by  Table  1.  Also  Figures  3,  5,  and  7 
shows  that  the  lower  the  stiffness  modification,  the  difference  between  the  time- 
steps  are  greater.  Consequently,  the  number  of  iterations  required  to 
approximate  the  response  is  lower  than  that  for  the  higher  stiffness  modifications. 
It  is  important  to  keep  in  mind  that  accuracy  is  not  just  determined  by  the 
stiffness  modification,  but  is  determined  by  the  time-step  as  well.  A  higher 
stiffness  modification  requires  a  smaller  time-step  for  accuracy  is  dependent  on 
the  ratio  of  the  stiffness  modification  to  the  time-step. 

It  should  be  noted  from  Table  1  with  a  0%  stiffness  modification  the  error 
is  very  low,  and  the  time-step  was  doubled  every  iteration.  Though  the  system 
was  approximated  in  seven  time-steps,  the  result  is  not  a  good  approximation  of 
the  system  response.  Though  Gordis  and  Neta  [3]  place  a  lower  limit  on  the  time- 
step  as  a  termination  requirement  in  case  a  function  cannot  be  approximated,  for 
systems  with  low  modifications  an  upper  limit  is  needed  in  order  to  accurately 
represent  the  system  response. 


22 


Finally,  the  measure  of  efficiency  needs  to  be  addressed.  The  ideal 
measure  of  efficiency  is  the  combination  of  speed  and  accuracy.  Clearly,  speed 
and  the  number  of  time-steps  are  directly  proportional  as  evidenced  by  Table  1. 
A  non-adaptive  trapezoidal  method  from  zero  to  two  seconds  at  a  time-step  of 
0.01  seconds  would  require  200  iterations  to  approximate  the  system  response. 
However,  accuracy  may  be  lost  when  using  a  non-adaptive  trapezoidal  method. 
Table  1  shows  that  for  a  stiffness  modification  of  10%  or  higher,  the  number  of 
iterations  is  greater  than  that  required  from  the  non-adaptive  approach.  Efficiency 
is  clearly  lost  in  terms  of  speed  with  an  adaptive  approach,  but  a  higher  degree  of 
accuracy  may  be  preserved. 

There  are  many  variables  that  have  to  be  taken  into  account  when 
discussing  efficiency.  Table  1  gives  the  results  for  the  most  obvious  variable  that 
can  affect  efficiency,  stiffness  modification.  Other  variables  to  consider  may  be, 
but  not  limited  to,  maximum  error  between  the  coarse  and  fine  approximations, 
and  the  acceptable  error  between  the  approximation  of  the  modified 
displacement  and  known  or  measured  data  (for  this  discussion  the  exact 
solution).  For  this  particular  model  it  is  shown  that  the  adaptive  approach  can 
give  reasonably  accurate  results  over  a  range  of  stiffness  modifications.  It  is  up 
to  the  analyst  to  consider  the  usefulness  of  this  method  for  their  particular  model. 

C.  SDOF  MASS-SPRING  SYSTEM  WITH  EXTERNALLY  APPLIED 

PERIODIC  FUNCTION  SUBJECTED  TO  A  STIFFNESS  MODIFICATION 

The  physical  representation  of  this  model  is  given  by  Figure  1  where  the 
external  excitation  is  now  a  sinusoidal  forcing  function.  The  stiffness  modification 
for  this  system  is  varied  in  the  following  intervals:  0%,  10%,  25%,  50%,  75%,  and 
100%  at  an  initial  time-step  of  0.01  second,  and  a  maximum  error  between 
coarse  and  fine  approximations  of  1x10  4 .  The  results  of  these  modifications  are 
presented  in  Table  2.  The  plots  of  modified  displacement  and  step  size  for  10%, 
50%  and  100%  have  been  included  to  show  the  variance  in  system  response  as 
the  stiffness  modification  is  raised.  Many  of  the  same  observations  for  this 
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system  follow  the  observations  made  in  the  previous  section  with  the  system 
subjected  to  an  external  step  excitation. 


Stiffness 

Greatest 

Least 

Total  number 

Time  to 

Modification(%) 

Error  (%) 

Error  (%) 

of  calculations 

complete(sec) 

0 

0.0 

0.0 

7 

0.1 

10 

0.8 

0.0 

638 

1.9 

25 

1.6 

0.0 

812 

2.8 

50 

3.1 

0.0 

1119 

3.9 

75 

3.8 

0.0 

1277 

4.5 

100 

4.6 

0.0 

1339 

4.7 

Table  2.  Response  to  stiffness  modifications  for  SDOF  mass-spring 
system  subjected  to  externally  applied  sinusoidal  excitation. 


An  interesting  point  to  note  is  that  the  error  obtained  from  this  system  is 
smaller  than  that  for  the  system  subjected  to  external  step  excitation  as  the 
stiffness  modification  is  increased.  This  result  would  make  sense  if  the 
corresponding  number  of  calculations  was  higher  for  the  periodic  function.  Since 
it  is  not,  the  reason  for  the  lower  error  is  not  understood  and  needs  to  be 
investigated  further. 


x*(t)  vs.  Time 


Figure  8.  Modified  response  for  a  SDOF  mass-spring  with  externally 
applied  sinusoidal  excitation  subjected  to  10%  stiffness 

modification. 
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Step  Size 


Step  Size  vs.  Step  Number 


Figure  9.  Step  size  for  a  SDOF  mass-spring  with  externally  applied  step 
excitation  subjected  to  10%  stiffness  modification. 


x*(t)  vs.  Time 


Figure  10.  Modified  response  for  a  SDOF  mass-spring  with  externally 
applied  sinusoidal  excitation  subjected  to  50%  stiffness 

modification. 
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Step  Size 


Step  Size  vs.  Step  Number 


Figure  1 1 .  Step  size  for  a  SDOF  mass-spring  with  externally  applied  step 
excitation  subjected  to  50%  stiffness  modification. 


x*(t)  vs.  Time 


Figure  12.  Modified  response  for  a  SDOF  mass-spring  with  externally 
applied  sinusoidal  excitation  subjected  to  100%  stiffness 

modification. 
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Step  Size  vs.  Step  Number 


Figure  1 3.  Step  size  for  a  SDOF  mass-spring  with  externally  applied  step 
excitation  subjected  to  100%  stiffness  modification. 


D.  MDOF  MASS-SPRING  SYSTEM  WITH  EXTERNALLY  APPLIED  STEP 
FUNCTION  SUBJECTED  TO  A  STIFFNESS  MODIFICATION 


The  physical  representation  of  this  model  is  given  by  Figure  14. 


Spring 

Modification 


Figure  14.  Two  DOF  mass-spring  system  subjected  to  stiffness 

modification. 


The  stiffness  modification  is  applied  to  the  second  spring,  and  the  external 
excitation  is  a  step  function  applied  to  the  second  mass  of  the  system.  Stiffness 
modifications  of  0%,  25%,  50%,  75%,  and  100%  are  made  at  an  initial  time-step 
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of  0.001  second,  and  a  maximum  error  of  MO  4  between  the  coarse  and  fine 
approximations.  The  results  can  be  found  in  Table  3. 


Stiffness 
Modification  (%) 

Greatest 
Error  (%) 
(Mode  1) 

Greatest 
Error  (%) 
(Mode  2) 

Total  number 
of  calculations 

Time  to 
complete(sec) 

0 

0.5 

0.5 

11 

0.08 

10 

1.7 

1.2 

3899 

91 

25 

2.1 

1.2 

5590 

181 

50 

1.7 

1.0 

7765 

343 

75 

2.0 

1.1 

8937 

450 

100 

1.8 

1.2 

9507 

505 

Table  3.  Response  to  stiffness  modifications  for  a  two  DOF  mass¬ 
spring  system  subjected  to  externally  applied  step  excitation. 


x*(t)  vs  'Exact' 


1  1 

A 

& 

\  \ 

♦  Modified  (Mode  1) 

♦  Modified  (Mode  2) 

- Unmodified  (Mode  1) 

- Unmodified  (Mode  2) 

- Exact  Modified  (Mode  1) 

- Exact  Modified  (Mode  2) 


0.8  1  1.2 


Figure  15.  Modified  response  for  a  two  DOF  mass-spring  with  externally 
applied  step  excitation  subjected  to  10%  stiffness  modification. 
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Figure  16.  Step  size  for  a  two  DOF  mass-spring  with  externally  applied 
step  excitation  subjected  to  10%  stiffness  modification 


x*(t)  vs  'Exact' 


Figure  17.  Modified  response  for  a  two  DOF  mass-spring  with  externally 
applied  step  excitation  subjected  to  10%  stiffness  modification. 


29 


Figure  18.  Step  size  for  a  two  DOF  mass-spring  with  externally  applied 
step  excitation  subjected  to  50%  stiffness  modification. 


x*(t)  vs  'Exact' 


Figure  19.  Modified  response  for  a  two  DOF  mass-spring  with  externally 
applied  step  excitation  subjected  to  100%  stiffness  modification. 
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Figure  20.  Step  size  for  a  two  DOF  mass-spring  with  externally  applied 
step  excitation  subjected  to  50%  stiffness  modification. 

From  Figures  15,  17,  19,  and  Table  3,  it  is  immediately  apparent  that  this 
approximation  has  a  very  low  error  and  remains  relatively  constant  as  the 
stiffness  modification  is  increased.  However,  the  number  of  calculations  per 
stiffness  modification  goes  up  as  well.  This  implies  a  smaller  time  step,  and 
therefore  accuracy  is  preserved  through  time-step  reduction.  This  can  be  seen  in 
Figures  16,  18,  and  20  where  the  time-steps  of  Figures  18  and  20  are  essentially 
half  the  value  of  the  time-steps  of  Figure  16. 

E.  MDOF  CANTILEVERED  ALUMINUM  BEAM  WITH  AN  EXTERNALLY 
APPLIED  STEP  EXCITATION,  SUBJECTED  TO  A  STIFFNESS 
MODIFICATION  TO  AN  ARBITRARY  BEAM  ELEMENT 

This  model  (Figure  21 )  is  the  most  important  model  of  this  work  because  it 
is  not  only  the  first  model  to  attempt  to  apply  the  adaptive  method  to  a  physical 
structure,  but  it  also  the  first  model  to  exhibit  conditional  stability.  As  mentioned  in 
the  Introduction,  the  magnitude  of  the  stiffness  modification  drives  the  integral 
equation  to  instability.  It  should  be  noted  that  stability  was  preserved  in  the 
previous  models  because  the  parameters  of  the  system  were  sufficiently  small. 
Once  conditional  stability  was  observed  in  this  model,  and  the  reason 
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determined,  the  parameters  of  the  SDOF  models  were  raised  to  the  magnitude  of 
the  cantilevered  beam.  The  programs  were  not  able  to  approximate  a  solution 
and  terminated  due  to  meeting  the  minimum  time-step  termination  requirement. 

Initially  the  goal  of  this  paper  was  to  apply  the  adaptive  method  to  the 
transient  analysis  of  a  structure.  Once  conditional  stability  was  exhibited  using 
the  non-adaptive  coarse  approximation,  the  goal  then  became  to  find  a  means  to 
restore  stability  to  the  approximation.  As  a  result,  the  half  solution  and  fine 
approximation  algorithms  were  not  developed. 

Figures  that  show  the  instability  for  stiffness  modifications  of  14%  and 
50%  are  provided  along  with  figures  in  which  appropriate  time-steps  have  been 
chosen  such  that  accurate  approximations  have  been  obtained.  The  14% 
stiffness  modification  was  chosen  because  below  this  value  the  instability  is  not 
obvious. 


Figure  21 .  Aluminum  cantilevered  beam  of  five  elements  with  an 
elemental  stiffness  change  to  the  third  element  subjected  to  a 
step  external  excitation. 
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x*(t)  (modified  displacement) 


x*(t)  vs.  Time 


Figure  22.  Instability  of  the  aluminum  cantilever  beam  with  a  stiffness 
modification  of  14%  and  a  time-step  of  0.01  second. 


Figure  23.  Instability  of  the  aluminum  cantilever  beam  with  a  stiffness 
modification  of  50%  and  a  time-step  of  0.01  second. 
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Figure  21  and  Figure  22  show  that  the  magnitude  of  the  instability  is 
directly  related  to  the  magnitude  of  the  stiffness  modification.  From  Equation  35  it 
is  obvious  that  as  the  stiffness  modification  is  increased,  the  magnitude  of  the 
kernel  increases,  and  therefore  the  approximation  becomes  more  unstable. 
Figure  23  and  Figure  24  show  that  a  step-size  of  sufficient  magnitude  can  restore 
stability  to  the  approximation.  The  key  observation  is  that  the  higher  the 
magnitude  of  the  stiffness  modification,  a  smaller  time-step  is  required  to  restore 
stability. 


x*(Q  vs.  Time 


Figure  24.  Stability  of  the  aluminum  cantilever  beam  with  a  stiffness 
modification  of  14%  and  a  time-step  of  0.0001  second. 
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x*(t)  vs.  Time 


Figure  25.  Stability  of  the  aluminum  cantilever  beam  with  a  stiffness 
modification  of  50%  and  a  time-step  of  0.00001  second. 

F.  THE  DERIVATIVE  OF  THE  TRANSIENT  STRUCTURAL  SYNTHESIS 
EQUATION  FOR  THE  MDOF  CANTILEVERED  ALUMINUM  BEAM  WITH 
AN  EXTERNALLY  APPLIED  STEP  EXCITATION,  SUBJECTED  TO  A 
STIFFNESS  MODIFICATION  TO  AN  ARBITRARY  BEAM  ELEMENT 

The  motivation  for  analyzing  the  derivative  of  the  synthesis  equation  is  to 
determine  if  the  derivative  is  unconditionally  stable.  Though  the  derivative  results 
in  the  velocity  of  the  system  response,  if  unconditionally  stable,  a  numerical 
integration  method  can  be  utilized  to  determine  the  position.  As  previously 
mentioned,  stability  of  integral  equations  is  a  complex  study  in  itself.  In  order  to 
test  stability  of  the  derivative  of  the  synthesis  equation,  the  same  time-step  that 
resulted  in  instability  were  used.  A  100%  stiffness  modification  was  implemented 
to  ensure  the  kernel  was  of  sufficient  size  to  induce  instability  in  the  original 
approximation  algorithm,  and  the  time  interval  was  over  the  time  that  the  system 
would  achieve  equilibrium.  The  system  was  also  run  with  the  same  parameters 
as  the  non-derivative  approximation  for  comparison  purposes. 
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From  Figure  26,  it  is  seen  that  at  100%  stiffness  modification,  0.01  second 
time-step,  the  system  exhibited  stability.  It  should  be  noted,  that  the 
approximation  began  to  trail  off  over  time.  Recall  in  the  Theory  section  that 
explained  in  detail  the  derivative  approximation,  the  modified  position  was 
derived  using  the  trapezoidal  rule  applied  to  velocities  that  were  derived  using 
the  trapezoidal  rule.  At  this  time-step,  the  error  is  large  and  can  accumulate 
quickly.  Any  accuracy  the  solution  may  have  had  is  lost  quickly. 


Non-adaptive  MDOF  (Derivative  Approach)(x'(t)) 


Figure  26.  MDOF  derivative  approximation  with  a  100%  stiffness 
modification  and  a  0.01  second  time-step  exhibiting  stability. 


Non-adaptive  MDOF  (Derivative  Approach)(x*(t)) 


Figure  27.  MDOF  derivative  approximation  14%  elemental  stiffness 
modification  and  a  time-step  of  0.01  second. 
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Non-adaptive  MDOF  (Derivative  Approach)(xT(t)) 


Figure  28.  DOF  derivative  approximation  50%  elemental  stiffness 
modification  and  a  time-step  of  0.01  second. 


It  can  be  seen  from  Figures  27  and  28  that  at  a  time-step  of  0.01  second, 
the  approximation  is  not  very  accurate,  but  it  exhibits  stability.  It  can  also  be  seen 
that  the  14%  stiffness  modification  most  closely  follows  the  exact  solution.  Now 
that  stability  has  been  achieved  using  the  derivative  approach,  the  next  question 
is  whether  the  time-step  can  be  refined  to  achieve  accuracy. 


Non-adaptive  MDOF  (Derivative  Approach)(x*(t)) 


Figure  29.  DOF  derivative  approximation  14%  elemental  stiffness 
modification  and  a  time-step  of  0.00001  second. 
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Non-adaptive  MDOF  (Derivative  Approach)(x*(t)) 


Figure  30.  DOF  derivative  approximation  50%  elemental  stiffness 
modification  and  a  time-step  of  0.00001  second. 


From  Figures  29  and  30,  it  can  be  seen  that  not  only  has  the  derivative 
approximation  achieved  stability,  but  it  also  accurately  approximates  the  system 
response.  For  the  stiffness  modification  of  14%,  step-size  of  0.00001  second,  the 
error  is  0.1%.  For  the  stiffness  modification  of  50%,  step-size  of  0.00001  second, 
the  error  is  0.7%.  Though  there  are  many  benefits  to  using  this  approach  to 
approximating  the  system  response,  the  time  needed  to  perform  the 
approximation  is  prohibitive  when  the  time-step  is  small.  This  is  the  reason 
Figures  29  and  30  are  plotted  over  a  small  time  interval. 

G.  THE  MATRIX  FORMULATION  FOR  THE  DERIVATIVE  OF  THE 
TRANSIENT  STRUCTURAL  SYNTHESIS  EQUATION  FOR  THE  MDOF 
CANTILEVERED  ALUMINUM  BEAM  WITH  AN  EXTERNALLY  APPLIED 
STEP  EXCITATION,  SUBJECTED  TO  A  STIFFNESS  MODIFICATION 
TO  AN  ARBITRARY  BEAM  ELEMENT 

The  final  point  of  the  previous  section  was  that  the  time  to  perform  the 
approximation  using  the  derivative  approach  was  not  desirable.  This  was  also 
true  for  the  non-derivative  approach  as  well.  Figures  25,  29,  and  30,  took 
approximately  two  hours  to  execute,  and  over  two  days  on  a  time-scale  from  zero 
to  two  seconds  in  order  to  capture  equilibrium  of  the  system. 
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The  time  that  it  takes  to  approximate  the  MDOF  systems  can  be  attributed 
to  the  algorithms  used  to  perform  these  approximations.  When  programming  the 
SDOF  systems  in  MATLAB,  there  are  many  time-saving  methods  that  avoid 
using  expensive  time  consuming  loops  when  operating  with  scalar  vectors.  When 
programming  these  algorithms  for  the  MDOF  cases,  the  same  time  saving 
methods  used  for  scalar  vectors  did  not  seem  to  have  a  direct  correlation  to 
matrix  and  vector  arrays.  As  a  result,  these  models  required  the  use  of  loops  to 
solve  the  approximations,  and  therefore  more  time  to  solve. 

The  matrix  formulation  for  the  derivative  approach  follows  the  formulation 
given  in  [1 1],  and  is  provided  in  Appendix  D.  The  elements  of  the  resulting  matrix 
are  matrices,  and  the  vectors  of  the  formulation  are  composed  of  vectors.  The 
challenge  with  the  matrix  formulation  is  that  in  order  to  approximate  the  solution 
on  a  time-step  needed  for  sufficient  accuracy,  the  matrix  is  large  and  requires  a 
large  amount  of  computer  memory  (RAM).  In  order  to  obtain  an  approximation, 
the  time  interval  of  approximation  needs  to  be  small.  There  are  supercomputers 
that  can  be  used  with  virtually  unlimited  memory  and  can  quickly  evaluate  large 
systems;  the  goal  is  to  develop  an  algorithm  that  can  be  evaluated  on  computers 
with  limited  memory,  a  common  restriction. 

The  challenge  is  to  determine  if  the  matrix  formulation  of  the  derivative  of 
the  synthesis  equation  can  produce  an  approximation  faster  than  the 
approximation  algorithms  previously  used.  A  comparison  of  the  matrix 
formulation  against  the  approximation  algorithms  over  small  time  intervals  was 
conducted  and  the  results  are  given  in  Table  4.  Table  4  gives  the  time  it  takes  to 
construct  the  components  of  the  matrix  formulation,  the  time  to  solve  the  system, 
and  the  time  it  took  to  execute  the  derivative  algorithm. 

From  Table  4  there  are  a  few  points  that  have  to  be  considered  prior  to 

determining  whether  the  matrix  formulation  is  a  viable  option.  First,  constructing 

each  part  of  the  matrix  formulation  takes  a  comparable  amount  of  time  to  the 

time  it  actually  takes  to  solve  the  system.  Constructing  the  matrix  formulation 

may  not  have  been  done  in  an  efficient  manner,  so  the  comparison  should  be 

39 


taken  from  the  point  at  which  the  fully  constructed  matrix  system  is  being  solved. 
Second,  on  small  time-intervals,  the  original  approximation  algorithms  are  faster 
than  constructing  and  solving  the  matrix  formulation.  It  is  difficult  to  tell  at  this 
point  as  to  whether  the  matrix  formulation  is  more  beneficial  for  time.  Due  to 
limited  computer  memory,  constructing  a  matrix  needed  to  solve  the  time  interval 
that  would  provide  a  good  comparison  could  not  be  accomplished.  Finally,  this 
method  should  not  be  fully  discounted.  The  results  that  it  produces  are  very 
accurate  (0.3%  error)  as  seen  in  Figure  31.  If  it  can  be  shown  that  on  a  longer 
time  scale  this  method  is  comparable  to  the  original  derivative  algorithms,  the 
matrix  formulation  can  be  a  valid  alternative  to  the  derivative  algorithms. 


Time 

Interval 

(sec) 

Construct 

Matrix 

Elements 

(sec) 

Construct 

Matrix 

(sec) 

Construct 

Vector 

(sec) 

Solving 

Matrix 

System 

(sec) 

Total 
time  to 
solve 
(sec) 

MDOF 

Derivative 

algorithm 

(sec) 

0-0.025 

20.1 

50.6 

0.12 

9.6 

83.1 

20.81 

0-0.050 

81.3 

210.0 

0.24 

67.8 

362.3 

77.3 

0-0.075 

183.1 

477.9 

0.35 

223.7 

883.7 

170.3 

Table  4.  Time  required  to  complete  each  part  of  the  matrix 
formulation  of  the  derivative  to  the  integral  formulation  for 
structural  synthesis. 


It  can  be  seen  from  Table  4  that  if  the  entire  process  of  constructing  and 
solving  the  matrix  formulation  is  the  measure  of  time  to  be  compared  against  the 
original  approximation  integral,  then  this  method  is  not  a  viable  alternative.  It 
should  be  pointed  out  that  the  total  time  for  the  matrix  formulation  given  in  Table 
4  is  an  average  of  three  runs.  These  runs  were  performed  to  check  for 
consistency.  However,  if  the  measure  of  time  is  the  time  that  it  takes  to  solve  the 
fully  constructed  matrix  system,  then  the  results  are  inconclusive.  Though  the 
matrix  formulation  is  faster  for  the  first  two  time  intervals,  on  the  last  interval  it 
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becomes  noticeably  slower.  Further  time  intervals  could  not  be  used  for 
comparison  due  to  the  limits  of  computer  memory. 

For  the  time  scale  of  0-0.075  seconds  with  a  time-step  of  MO  5  sec,  the 
matrix  size  is  30000x30000,  a  relatively  large  system  to  be  solved.  Though  the 
program  was  able  to  approximate  the  response  for  a  time  interval  of  0  to  0.10 
seconds,  computer  performance  was  severely  affected  and  the  time  data 
gathered  was  too  varied  to  establish  a  trend  for  comparison.  Though  the  time 
data  was  too  varied  to  be  used,  Figure  31  gives  the  result  for  the  matrix 
formulation.  It  can  be  seen  that  this  method  is  incredibly  accurate  with  an  error  of 
0.3%. 


Matrix  Formulation  of  the  Derivative  of  the  Synthesis  Equation 


Figure  31 .  Matrix  formulation  of  the  derivative  of  the  synthesis  equation 
for  the  Aluminum  cantilevered  beam  with  a  100%  stiffness 
modification  and  time-step  of  .00001  second. 
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IV.  CONCLUSION  AND  RECOMMENDATIONS 


The  adaptive  method  to  solving  the  integral  formulation  for  transient 
structural  synthesis  is  a  valid  approach.  As  of  now,  the  adaptive  method  is  best 
used  to  approximate  the  transient  response  of  SDOF  and  two  DOF  systems.  In 
terms  of  the  SDOF  system,  the  programming  methods  available  bypass  the  need 
to  use  time  consuming  loops  allowing  the  approximations  to  be  performed 
quickly.  In  terms  of  the  two  DOF  system,  this  particular  MDOF  system  is 
sufficiently  simple  enough  that  even  though  the  use  of  loops  is  necessary,  the 
approximations  are  performed  with  sufficient  speed. 

Due  to  computational  inefficiency,  the  adaptive  method  is  not  an  efficient 
alternative  to  solving  MDOF  transient  structural  synthesis  problems  over  a  long 
time  scale.  The  next  step  is  to  study  the  programming  methods  available  that  can 
efficiently  perform  the  adaptive  algorithms  using  arrays  of  vector  and  matrices. 
Once  these  methods  are  developed,  then  adaptivity  can  be  implemented. 

A  matrix  formulation  of  the  MDOF  problem  was  considered  to  overcome 
the  computational  inefficiency  due  to  using  array  operations.  Unfortunately,  the 
size  of  matrix  needed  to  approximate  the  response  of  the  system  is  very  large 
and  limited  by  available  computer  memory.  Approximations  can  only  be 
performed  on  small  time  intervals  due  to  available  memory,  but  on  these  small 
intervals,  the  adaptive  approximation  algorithms  are  comparable  in  speed.  Since 
the  analysis  is  assumed  to  be  conducted  by  computers  with  limited  memory,  an 
efficient  method  such  as  a  block-by-block  method  needs  to  be  investigated.  A 
conclusive  comparison  of  the  matrix  formulation  to  the  adaptive  approximation 
algorithms  cannot  be  made  at  this  time. 

If  the  structural  modification  is  sufficiently  large,  then  the  integral 
formulation  for  transient  structural  analysis  becomes  conditionally  stable.  The 
derivative  of  the  integral  formulation  for  transient  structural  synthesis  restores 
unconditional  stability  and  provides  velocity  of  the  modified  system  response. 
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Keep  in  mind  that  the  derivative  of  the  synthesis  equation  complicates  the 
adaptive  algorithms. 

Finally,  once  computational  efficiency  for  MDOF  problems  has  been 
established,  mass  and  damping  modifications  need  to  be  studied.  Once  these 
modifications  are  understood  and  successfully  approximated,  then  consideration 
needs  to  be  given  to  non-linear  structural  modifications. 


44 


APPENDIX  A.  COARSE  SOLUTION  FOR  THE  INTEGRAL 
EQUATION  FORMULATION  FOR  STRUCTURAL  SYNTHESIS 


The  MDOF  integral  formulation  of  structural  synthesis  is  given  by  the  following 
equation: 


i 

jx*(f)}  =  {x(t)j  -  J [h(t  -  r)][A^]|x*(r)| dr 


(43) 


At  time-step  / ,  the  integral  can  be  written  as  the  sum  of  integrals: 


-1  <M 


x  }  =  J  -r)][AS']|x*|  dr  + 


(44) 


Applying  the  trapezoidal  rule: 

J  i= 0  z 


X  }  + 

+1 


[/iL-[AK']{x*})  +  {x}. 


(45) 


}, = [h]hM  [hlj  MK1.+H  (46) 

J  i= 0  z  i= 0  z 


% = ~i\si  [hl,i  Wjj  (47) 

J  i= 1  z  i= 0  z 


m,=-£^[*umm, 

J  1=0  z 

Zj  [A],,  M  {V},  -  ^  [h]..  [AK]  { x}  +  i  x 

1=0  z 


2  1  L 


y  1  ^ 


(48) 


H.,  MK),=T^  MM, +{*'}, 


i= 0 


(49) 


a([^]  +  ^;[ALM){*T=_|i(‘S'+^‘)[',LMK}, 


+  {x 


(50) 


From  Equation  (1),  it  can  be  seen  that  when?  =  0,  |x*}q  =  {x}0  ={0} . 
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APPENDIX  B.  THE  DERIVATIVE  OF  THE  INTEGRAL  EQUATION 
FORMULATION  FOR  STRUCTURAL  SYNTHESIS 


Starting  with  the  MDOF  integral  formulation  of  structural  synthesis: 

t 

jx* (r)j  =  |x(7)}  -  J [h(t  ~  r)] [ ] jx*(r) j  d v 


Taking  the  derivative  of  Equation  (1)  with  respect  to  t: 

d  (  * ,  A )  d  (  ,  A  ^  d  f 

dt 


*  ^ = ^  ^  i  [/?(r  ~  r)]  iak  d  z 


Recall: 


{-<t)}  =  \[h(t-r)]{f(T)}di 


Substituting  Equation  (3)  into  Equation  (2): 


d_ 

dt 


x  (0}  =tj\ [Kt  -  r)] {f(j)}d t  - y  J [h(t  -  r)] [ AK] {x* (r)}  d r 


d  r, 


dt 


dt  • 


Analyzing  —  {x(r)}  : 


d  r 


{*(0}  =  —  J  [h(d  ~  r)]{f(r)}d  i 


d  r 


\[h(t  -  T)]{f(T)}dT  =  [/!((' -  0]  {/(0} ^ - [h(t -0)]{/(0)} 


dt 


dt 


+JT([/?(r_r)]{^(r)})Jr 


d(  0) 
dt 


d  r 


dt 


J  [hit  -  r)]  {f(r)}d  r  =  0  -  0  +  J  [/?(>  -  r)]  {fir)}  d  r 


d  r 


dt 


J  [Kt  ~  r)]  {fir)}d r  =  {x(r)}  =  J  [h{t  -  r)]  {/(r)}  dr 


d  r. 


Similarly  when  analyzing  —  J[/i(r-r)][M']{x*(r)j<ir : 


(51) 

(52) 

(53) 

(54) 

(55) 

(56) 

(57) 

(58) 
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(59) 


d_ 

dt 


i  i 

J[/7(f-r)][A/T]  jx*(r)jrfr  =  J  h(t  -  z)  [A/T]{x*(r)j  dr 


Substituting  Equation  (8)  and  Equation  (9)  into  Equation  (2)  gives: 

t  t 

|x*(r)|  =  J[/i(r-r)  {/(r)}<ir- J[^i(r-r)  [A^]jx*(r)j dz 


(60) 


which  reduces  to: 


i 

x*(?)}  =  {x(r)}-J|^/i(r-r)  [A^r]|x*(r)| dr 


(61) 


It  can  be  seen  that  when  r  =  0,  |x*(0)|  ={x(0)}  =  {0}.  Also  recall  that  from 
Equation  (1),  |x*(0)|  =  {x(0)}  =  {0}  when  r  =  0. 
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APPENDIX  C.  COARSE  SOLUTION  TO  THE  DERIVATIVE  OF  THE 
INTEGRAL  EQUATION  FORMULATION  FOR  STRUCTURAL 

SYNTHESIS 


Starting  with  the  derivative  of  the  integral  formulation  for  structural  synthesis: 


i 

x*(r)}  =  {i(r)}-||^/i(r-r)  [A^]jx’f(r)j dz 


The  integral  is  broken  into  the  sum  of  integrals; 


{x*(r)}  =  J  -  r)][AAT]  {**(*■)}  dr  +  {i(f)} 


Applying  the  trapezoidal  rule: 

,  .  iz1  i 


(=0  z 


J  i=0  z  ■ 


{iT=-Zv5«(['iL[AX]{/},)-Z^.,(W,,M{z}> 

i=0  Z  i= 0  z 


The  unknown  |x*|  needs  to  be  expressed  in  terms  of  the  derivative  terms  by 


means  of  the  trapezoidal  rule: 

,  ,  1 


Substituting  Equation  (7)  into  Equation  (6): 


J^](Z^(K}i+{i'L)  +  7<J,(K},,+{i'},)  +  {Z}0)  +  WJ 
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,  jM  i 

-t‘  =-Y- 


-It(<5,+^,)[a].1[aa:]{x-} 

i=0  z  y’ 

\s>  [\L  ^ 


-kM„M<£^4*VKL» 


-JSJ  M  „  [“]{*'};.,  ~\S>  [*],,  MM,  +  {*}, 
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APPENDIX  D.  THE  MATRIX  FORMULATION  TO  THE  DERIVATIVE 
OF  THE  TRANSIENT  STRUCTURAL  SYNTHESIS  EQUATION 


In  order  to  derive  the  matrix  formulation,  it  is  helpful  to  look  at  a  number  of 
iterations.  The  derivative  of  the  transient  structural  synthesis  equation  takes  the 
form: 


Let  t=0: 


Lett  =1: 


{*•}-{ 

k. 

TT 

* 

< , 
i  i 

TT 

i 

^  | 

%~T' 

(73) 

{*«>)}  = 

0 

|i(0)}- J^/i(0)J[A^]|x*(r)|  dr 

0 

(74) 

•'■{^*(0)}  ={i(0)} 

(75) 

{*'«}  =  { 

1 

x*  (1)  J  -  J  \ji(\  -  t)  J  [AK  ]  |  jc*(r)  Jrfr 

0 

(76) 

(i(l)}-Ar 

}WJAJfK/(0))4Wn[AJfIl/(1>!) 

(77) 

{.f(l)}  =  {.v(l)}- |-[/i]]o  [AK]{*(0)} ~[h\  [AAT]f  ^{i'(0)}  +  ^{**(1)}  +  {x*(0)} )  (78) 


t(l)}  = 


2  L  Jio 
A  t 


2  L  -'ll 1 

A]11[AA’]!i’<o>} 


(79) 


1  + 


At2 


a]„[a^]  j.i-(D} 


Following  the  same  process  for  t=2,  and  t=3  gives: 
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{.i(2)}  =  (  ^[a]m  [AK-]+ to[h\  [AK]+  y  [A]2;  [A  K]  +  l{*(0)} 


At 

~2 


f  At2 


f  At2 


'■l.M+y-lALM  {*(»>} 


v  ^ 
f 

1  + 


AT 

4  L  J22  1 

/ 


At 


/,]22[AX]  [i*(2)} 


4  L  J22 


{*(3)}  = 


(At2 


>]„[**]■  2 


+ 


A'  f/,L  [A/f]  +  ^-[*]33  [AAT]){a-*(°)} 

A 


^[*]„  M  +  A»2  [h\  [A  K]  +  ^\  [AAT]j{i'(l)} 

A 


V 

t[*]»[ak]+ 


At 

~2 


-IaLmV®} 

J 


'•lm'ik®} 


Let  t  =i: 


{*(;')}  = 


(80) 


yWj.[Ai:]+A»[,*],1  [ak]+a»[a]m  [AAr]+y[*]3j  [Ai:]+j{.r'(0)} 

'  ,  .9  .  ?  .2  'A 


[*L  M + ^  ['■y  + - + y  Mi-  M +)  { /(0)) 

r^M2l[A«]+y1['-],2[A«]+--+^['-]i[A«]){-*'(l))} 


‘';]J1[A«]+y1W/JA«]+- 


V 

"at2- 

l  2 


L-Jji 


L  Ji? 
At 


lU'i  1 


(81) 


At2  ^ 

[A K]  +  At2  [h\  ,2  [A K]  +  ...  +  —  [/z]  [A^r]  {i*(l)}  (82) 


At  r  • 


^\Aak\  {*(^} 


1  + 

V  4  L  Jjj  j 


Assuming  rest  initial  conditions: 

WO)}  =  {x(0)}  =  {x*(0)}  =  {i*(0)}  =  {0} 


Putting  the  system  of  equations  into  matrix  form: 


52 


[0] 


[0] 


[0] 


Equation  83  is  the  equation  to  be  solved. 
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APPENDIX  E.  MATLAB  CODE  FOR  A  SDOF  MASS-SPRING 
SYSTEM  WITH  EXTERNALLY  APPLIED  STEP  FUNCTION 
SUBJECTED  TO  A  STIFFNESS  MODIFICATION  (SUPPORTING 

FUNCTIONS  INCLUDED) 


S-S-S-S-S-S-S-S-S'S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  is  the  main  program  for  an  adaptive  SDOF  mass-spring  problem  % 
%subjected  to  a  periodic  forcing  function.  % 

%  Created  by  Keenan  L.  Coleman  8/25/2014  % 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


close  all; 
clear  all; 
clc; 

%  Intializing  the  SDOF  parameter  for  the  exact  solution: 

damp  ratio  =0.1;  %  damping  ratio 

k  =  100;  %  spring  constant  (N/m) 

delk  =  .  l*k;  %  stiffness  modification 

mass  =0.1;  %  mass  (kg) 

omega2  =  sqrt((k  +  delk) /mass);  %  natural  frequency  due  to 

stiffness  modification  (used  to  calculate  the  exact  solution) 
dampfreq2  =  omega2*sqrt ( 1-damp  ratioA2);  %  damped  frequency 
stiffness  modification  (used  to  calculation  the  exact  solution) 
%  Initializing  the  adaptive  algorithm  variables: 

%  Time  variables: 

%  start  time  (sec) 

%  end  time  (sec) 

%  initial  time-step/time  difference 


(sec) 

minimum  time  step  that  terminates  the  adaptive 


time_start  =  0; 
time  end  =  2; 
de 1 t  =  .01; 
min  del  =  0.00001; 
process  during  step  halving 

t  =  time  start :min  del: time  end;  %  time  interval  for  the  function 
handle  to  the  exact  solution  to  the  modified  response 
t  step  =  [0  delt] ;  %  initializing  the  vector  of  time-steps 

t  diff  =  [delt  0];  %  initializing  the  vector  of  time  differences 

between  time-steps  for  the  kernel 

t  stepcount  =  time  start  +  delt;  %  the  value  of  the  current  time- 
step 


h  =  [0  (delt  -  time  start)];  %  initializing  the  vector  of  time 

differences  for  the 

h2  =  [0  (delt  -  time_start) ] ; 

%  Spacial  variables: 
x  initial  =  0; 
x_coarse  =  [x_initial] ; 
x  half  =  [x  initial]; 
x  fine  =  [x_initial]; 
x  approx  =  [x  initial]; 
count  =  1; 
count2  =  0; 
double_count  =  0; 
half^count  =  0; 

max  error  =  .00001;  %  maximum  error  between  fine  and  coarse  solution 

to  determine  step-halving  or  step-doubling 
%  The  exact  solution... 
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exact_solution  =  @ (t) (1/ (mass*dampfreq2) ) *exp (- 

damp  ratio*omega2*t) . * (dampf req2*exp (damp  ratio*omega2*t) - 

damp  ratio*omega2*sin (dampfreq2*t) - 

dampfreq2*cos (dampfreq2*t) ) / (damp  ratioA2*omega2A2+dampfreq2A2) ; 

exact_solution  =  exact_solution (t) ; 

tic 

%  Solving  for  x(t),  the  unmodified  solution  (uses  the  Trap  rule 
solution) . 

[x  solution]  =  x  solve_2 (time  start,  time  end,delt,  min  del); 
while  t_stepcount<=time  end 

%  Calling  values  for  x(t) (  This  is  g(t)  in  Dr.  Neta's  paper) . 

[x  result,  x  result  half]  =  x  solution  call  2 (time_start, 
min_del, time_end, t_step (count)  , t_step (count+1) , x_solution) ; 

%  This  is  the  kernel  function 

[  kern  elem, kernl , kernhalf  ]  =  G2  kernel (  t_diff , t_step, count, delt 

)  ; 

%  Course  Solution 
[x_coarse,  x^coarsevalue]  = 

G  course_solution^2 (h, x  coarse,  x  result,  kern  elem, kernl , count) ; 

%  x-half  solution 
[x  half,  x  halfvalue]  = 

G  half_solution  2 (h, x^coarse, x  half , kernhalf , kernl , count, x_result  half) 

r 

%  Fine  Solution 
[x  fine,  x  finevalue]  = 

G  fine  solution  2 (h, x^coarse, x  fine,x  halfvalue, x_result, kern  elem, kern 
1 , kernhalf,  count)  ; 

%  evaluating  the  error 

x  value  =  [x  finevalue, x  coarsevalue] ; 

relerror  =  abs (dif f (x^value) ) /abs (max (x_coarse) ) ; 

relerror  =  0 . 99*relerror; 

inti  =  f ind (t==t_stepcount)  ; 

t_stepcount; 

exact  =  exact_solution ( inti ) ; 

x  finevalue; 

x_coarsevalue; 

relerror; 

%  adjusting  the  time-step 
if  delt  <=  min  del 

sprintf (' Minimum  time  difference  exceeded') 
break 

end 

if  relerror  <  max^error 
x  approx (count+1 )  =  x  finevalue; 
x  coarse (count+1 )  =  x_finevalue; 
if  t  stepcount  ==  time  end 
break 

end 

count  =  count  +1; 
count2  =  count2  +  1; 
double_count  =  double_count  +  1; 
delt  =  2  *  de 1 1 ; 
h (count+1)  =  delt; 
h2(count2+l)  =  delt; 

t_step (count+1 )  =  t_step (count)  +  delt; 
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t_stepcount  =  t_step (count+1 ) ; 
t^stepvalue  =  t_step (count)  +  delt; 

[rl,cl]  =  size (t_step) ; 
for  ii  =  l:cl 

t_diff(ii)  =  t_step(end)  -  t_step(ii) 

end 

end 

if  relerror  >  max^error 
delt  =  de 1 t / 2 ; 
count2  =  count2  +  1; 
half_count  =  half_count  +  1; 
h (count+1)  =  delt; 
h2(count2+l)  =  delt; 

t_step (count+1 )  =  t_step (count+1 )  -  delt; 
t_stepcount  =  t_step (count+1 ) ; 
t_stepvalue  =  t_step (count+1 )  ; 

[rl,cl]  =  size (t_step) ; 
for  ii  =  l:cl 

t_diff(ii)  =  t_step(end)  -  t_step(ii) 

end 

x_coarse  =  x^coarse (1 : count) ; 
x  half  =  x  half (1 : count) ; 
x  fine  =  x  fine (1 : count) ; 

end 

if  t  stepcount  ==  time  end 
break 

end 


if  t  stepvalue  >  time^end 

t_step (count+1 )  =  time_end; 

h(count  +1)  =  time_end  -  t_step (count) ; 

h2 (count2  +1)  =  time_end  -  t  step (count); 

t_stepcount  =  t_step(count  +1); 

delt  =  h (count+1) ; 

[rl,cl]  =  size (t_step) ; 
for  ii  =  l:cl 

t_diff(ii)  =  t_step(end)  -  t_step(ii) 

end 

x^coarse  =  x^coarse (1 : count) ; 
x  half  =  x  half (1 : count) ; 
x  fine  =  x  fine (1 : count) ; 

end 


end 

toe 

double_count 
half  count 

total_count  =  count2 
%  Plotting  the  results: 

%  Change  default  axes  fonts, 
set ( 0 ,  ' Def aultAxesFontName ' , 
set(0,  ' Def aultAxesFontSize '  , 
%  Change  default  text  fonts. 
set(0,  ' Def aultTextFontname '  , 
set(0,  ' Def aultTextFontSize '  , 
figure  ( 1 ) 


' Arial ' ) 

14) 

' Arial ' ) 

14) 
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set (gcf , ' WindowStyle ' , 'Docked') 

plot (t_step, x_approx, ' kd  , t, exact_solution, ' k ' ) 
legend ( ' x* (t) ' , ' Exact ' ) 
xlabel('Time  (sec) ') 

ylabel ( ' x* (t)  (modified  displacement)') 
title('x*(t)  vs.  Time') 
grid  on 
figure (2 ) 

set (gcf , 'WindowStyle', 'Docked') 

[rl,cl]  =  size(h2); 

tt  =  l:cl; 

plot (tt, h2 , ' bd ' ) 

xlabel ( ' Step  Number') 

ylabel (' Step  Size') 

title (' Step  Size  vs.  Step  Number') 

grid  on 

axis ( [0, total_count,  0,  . 045] ) 
figure  ( 3 ) 

set (gcf , ! WindowStyle  ,  Docked') 

[rl,cl]  =  size(h2); 

tt  =  l:cl; 

plot (tt, h2 , ' bd ' ) 

xlabel (' Step  Number') 

ylabel (' Step  Size') 

title (' Step  Size  vs.  Step  Number') 

grid  on 

axis ( [0, 150,0, .045]) 

%  Evaluating  the  error... 

exact_solution_error  =  @ (t_step) (1/ (mass*dampfreq2) ) *exp (- 

damp  ratio*omega2*t  step) .* (dampfreq2 *exp (damp  ratio*omega2*t_step) - 

damp  ratio*omega2*sin (dampfreq2*t  step) - 

dampfreq2 . *cos (dampfreq2*t_step) ) / (damp_ratio A2 *omega2 A2+dampf req2 A2 ) ; 

exact_solution_error  =  exact_solution_error (t_step) ; 

max  diff  =  max (abs (exact  solution  error  -  x  approx)); 

min^diff  =min (abs (exact  solution  error  -  x  approx)); 

error  eval  =  max (abs (exact  solution  error  - 

x  approx) ) /max (exact_solution  error) 

error  eval2  =  min (abs (exact  solution  error  - 

x  approx) ) /max (exact  solution  error) 


S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  is  the  coarse  solution  algorithm  as  derived  by  Gordis  and  Neta. 
%Created  by  Keenan  Coleman  on  8/25/2014 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

function  [x^coarse, x^coarsevalue]  = 

G_course  solution^2 (h, x  coarse, x  result, kern  elem, kernl , count) 

lhsCoef f^coarse  =  1  -  0 . 5*h (end) *kernl (end) ; 
for  ss  =  1 : count 

h_term(ss)  =  (h(sstl)  +  h(ss)); 

end 

x_product  =  0 . 5*h_term. *kern  elem. *x_coarse; 
x  sum  =  sum (x_product)  +  x  result; 
xcoarse (count+1 )  =  x_sum/lhsCoef f_coarse; 
x  coarsevalue  =  x  sum/lhsCoeff  coarse; 
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%This  is  the  half  solution  algorithm  as  derived  by  Gordis  and  Neta% 
%Created  by  Keenan  Coleman  8/25/2014  % 
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ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

function  [x  half,x  halfvalue]  = 

G_half  solution  2 (h, x^coarse, x  half , kernhalf, kernl , count, x  result_half) 
lhsCoef f_half  =  1  -  0 . 25*h (end) * kernl (end) ; 
x_product  half  =  0; 

if  count  >=  2; 

for  ss  =  l:count-l 

h^term(ss)  =  (h(ss+l)  +  h(ss)); 

end 

x_product_half  =  0.5*h  term. *kernhalf ( 1 : count  - 
1 ). *x_coarse ( 1 : count  -1)  ; 
end 

x  sum  =  sum (x_product_half )  +  0 . 5* (h (count)  +  0.5*h (count  + 

1 )) *kernhalf (count) *x_coarse (count)  +  x  result  half; 
x  half (count+1 )  =  x  sum/lhsCoeff  half; 
x  halfvalue  =  x  sum/lhsCoeff  half; 

S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  is  the  half  solution  algorithm  as  derived  by  Gordis  and  Neta% 
%Created  by  Keenan  Coleman  8/25/2014  % 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


function  [x  half,x  halfvalue]  = 

G_half_solution_2 (h, x^coarse, x  half, kernhalf, kernl , count, x_result_half ) 
lhsCoef f_half  =  1  -  0 . 25*h (end) *kernl (end) ; 
x_product  half  =  0; 

if  count  >=  2; 

for  ss  =  l:count-l 

h_term(ss)  =  (h(ss+l)  +  h(ss)); 

end 

x_product  half  =  0.5*h  term. *kernhalf ( 1 : count  - 
1 ). *x_coarse ( 1 : count  -1); 
endd 

x  sum  =  sum (x_product  half)  +  0 . 5* (h (count)  +  0.5*h (count  + 

1) ) *kernhalf (count) *x_coarse (count)  +  x_result  half; 
x  half (count+1 )  =  x  sum/lhsCoeff  half; 
x  halfvalue  =  x  sum/lhsCoeff  half; 

S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%  This  calculates  the  kernel  at  the  time-steps  and  half  time-steps% 

%as  required  by  the  algorithms  derived  by  Gordis  and  Neta.  % 

%Created  by  Keenan  Coleman  on  8/25/2014  % 
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ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

function  [  kern  elem, kernl , kernhalf  ]  =  G2  kernel ( 
t^diff , t_step, count,  delt  ) 

%  intializing  the  SDOF  parameters... 
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damp  ratio  =  0.1; 
k  =  100; 
delk  =  .  1  *  k ; 
mass  =  0.1; 

omegal  =  sqrt ( k/mass ) ; 
a  =  damp  ratio*omegal ; 

dampfreq  =  omegal*sqrt(  1  -  damp  ratioA2); 

[rowl , columnl ]  =  size  (t  diff); 
for  ii  =  1: (columnl-1) 

t_half(ii)  =  (t_step (count) +t_step (count+1 )) /2  -  t_step(ii); 

end 

kern  =  @ (t  diff) (-delk  / (mass*dampfreq) ) *exp (- 

a*t_diff ) . *sin (dampf req*t_dif f ) ; 

kernl  =  kern(t  diff); 

kernhalf  =  kern(t  half ) ; 

kern^elem  =  kernl (1 : end-1) ; 

end 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%  This  function  calls  the  value  for  the  pre-calculated  baseline 
%response . 

%Created  by  Keenan  Coleman  on  8/25/2014 

S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

function  [x_result,x  result  half]  =  x  solution  call  2 (time_start, 
min  del, time  end,  t  stepl,  t  step2 , x^solution) 

xx  =  time  start :min  del: time  end; 

cl  =  find(xx  >=  t_step2); 

c2  =  find(xx  >=  (t_stepl  +t_step2 ) /2 ) ; 

cl  =  min (cl ) ; 

c2  =  min (c2 ) ; 

x_result  =  x^solution (cl ) ; 
x  result  half  =  x  solution (c2 ) ; 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'9'9'9'9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  function  calculates  the  baseline  response  whose% 

%values  will  be  called  as  needed.  % 

%Created  by  Keenan  Coleman  on  8/25/2014  % 

S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooo 

function  [x  solution]  =  x  solve_2 (time_start,  time  end,delt,  min  del) 

damp_ratio  =  0.1; 
k  =  100; 
delk  =  . 1 *  k ; 
mass  =  0.1; 

omegal  =  sqrt ( k/mass ) ; 
a  =  damp_ratio*omegal ; 

dampfreq  =  omegal*sqrt(  1  -  damp  ratioA2); 
countl  =  1; 
count2  =  2; 

gi  ( 1 )  =  0; 

t  =  time_start : delt : time  end; 

%  Implementing  the  algorithm: 

for  jj  =  time_start  +  delt : delt : time  end 

for  ii  =  l:count2 
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fl (ii)  =  1/ (mass*dampf req) *exp (-a* ( j j -t (ii) ) ) *sin (dampfreq* ( j j - 
t ( ii )  )  )  ; 
end 

fl(l)  =  . 5*f 1 (1) ; 
fl (end)  =  . 5*fl(end); 
g_value  =  delt*sum ( f 1 ) ; 
gl(countl+l)  =  g_value; 
countl  =  countl  +  1; 
count2  =  count2  +  1; 
end 

xx  =  time  start :min  del: time  end; 
yy  =  spline (t, gl, xx) ; 
x_solution  =  yy; 

exactunmod  =  1/ (mass*dampfreq) *exp (-a*t) . * (dampfreq*exp (a*t) - 
a*sin (dampfreq*t) -dampf req*cos (dampfreq*t) ) / (aA2  +  dampfreqA2); 
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APPENDIX  F.  MATLAB  CODE  FOR  A  SDOF  MASS-SPRING 
SYSTEM  WITH  EXTERNALLY  APPLIED  PERIODIC  FUNCTION 
SUBJECTED  TO  A  STIFFNESS  MODIFICATION  (SUPPORTING 

FUNCTIONS  INCLUDED). 


ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%  This  models  a  SDOF  problem  with  a  periodic  forcing  function. % 
%  Created  by  Keenan  Coleman  on  8/25/2014  % 

S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


close  all; 
clear  all; 
clc; 

%  Initializing  variables: 

delt  =  .01;  %  Initial  step  size 

t_start  =0;  %  Start  time 

t_stop  =2;  %  Stop  time 

time  =  t_start  :  .00001  :  t  stop  ;%  Timeline 
nStep  =  length (time) ;  %  Number  of  steps 

min_step  =  .00001;  %  Minimum  step  size  (necessary  for 

method) 

%  Physical  Parameters: 


m  =  0.1; 

k  =  100; 

wn  =  sqrt (k/m) ; 

fn  =  wn/ (2*pi)  ; 

z  =  0.01; 

wd  =  wn  *  sqrt(l 

delk  =  1  *  k; 


%  Mass 

%  Spring  Stiffness 
%  Natural  frequency  in  rads/sec 
%  Natural  frequency  in  Hz 
%  Damping  ratio 
zA2);%  Damped  frequency 

%  Change  in  stiffness 


adaptive 


%  Evaluating  the  baseline  response: 
tic 

[  x_base  ]  =  x_baseline (t__start,  t_stop,  min^step,  m,  wn,  wd,  z) ; 
%  Implementing  the  adaptive  algorithms  from  the  Neta  paper: 

%  Initializing  the  algorithm  variables: 

max  error  =  .0001; 

t_step  =  [0  delt] ; 

tdiff  =  [delt  0]; 

x  initial  =  0; 

h  =  [0  (delt  -  t_start) ] ; 

h2  =  [0  (delt  -  t_start) ] ; 

x_coarse  =  [x_initial] ; 

x  half  =  [x  initial]; 

x  fine  =  [x_initial]; 

x  approx  =  [x_initial]; 

count  =  1; 

t_stepcount  =  t_start  +  delt; 

count  =  1; 

count2  =  0; 

double_count  =  0; 

half  count  =  0; 
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tic 

while  t_stepcount<=t_stop 

%  Calling  the  x  base  and  half  x  base  values  needed  for  the 
coarse, half 

%  and  fine  values. 

[base_result, base_result_half ]  =  x_base_call (t_start, 
min_step, t_stop, t_step (count)  , t_step (count+1 ) , xbase) ; 

%  Evaluating  the  necessary  kernel  values: 

[  kern  elem, kernl , kernhalf  ]  =  kernel  eval (  m,wn,wd,z, 
t_diff , t_step,  count  ); 

%  Evaluating  the  coarse  solution: 

[x_coarse, x_coarsevalue]  = 

x_course_solution (h, x^coarse, base  result, kern^elem, kernl , count, delk) ; 

%  Evaluating  the  half-solution : 

[x  half,x  halfvalue]  = 

x  half  solution(h,x  coarse, x  half , kernhalf , kernl , count, base  result  half 

, de 1 k ) ; 

%  Evaluating  the  fine  solution: 

[x  fine,  x  finevalue]  = 

x  f ine__solution (h, x^coarse, x  fine,x  halfvalue, base  result, kern  elem, ker 
nl ,  kernhalf,  count,  delk)  ; 

%  evaluating  the  error 
if  x_coarse==0 

relerror  =  0; 

else 

x  value  =  [x  finevalue, x  coarsevalue] ; 

relerror  =  abs (dif f (x^value) ) /abs (max (x_coarse) ) ; 

relerror  =  0 . 99*relerror; 

end 

%  adjusting  the  time-step 
if  delt  <=  min  step 

sprintf (' Minimum  time  difference  exceeded') 
break 

end 

if  relerror  <  max^error 
x  approx (count+1 )  =  x  finevalue; 
if  t_stepcount  ==  t_stop 
break 

end 

count  =  count  +1; 
count2  =  count2  +  1; 
double_count  =  double_count  +  1; 
delt  =  2  *  de 1 1 ; 
h (count+1)  =  delt; 
h2(count2+l)  =  delt; 

t_step (count+1 )  =  t_step (count)  +  delt; 
t_stepcount  =  t_step (count+1 )  ; 
t_stepvalue  =  t_step (count)  +  delt; 

[rl,cl]  =  size (tstep) ; 
for  ii  =  l:cl 

t_diff(ii)  =  t_step(end)  -  t_step(ii); 

end 

end 
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if  relerror  >  iax_error 
de 1 1  =  del t / 2 ; 
count2  =  count2  +  1; 
half_count  =  half_count  +  1; 
h(count+l)  =  delt; 
h2(count2+l)  =  delt; 

t_step (count+1 )  =  t_step (count+1 )  -  delt; 
t_stepcount  =  t_step (count+1 )  ; 
t_stepvalue  =  t_step (count+1 ) ; 

[rl,cl]  =  size (tstep) ; 
for  ii  =  l:cl 

tdiff(ii)  =  t_step(end)  -  t_step(ii) 

end 

xcoarse  =  x_coarse (1 : count) ; 
x  half  =  x  half (1 : count) ; 
x  fine  =  x  f ine (1 : count) ; 

end 

if  t_stepcount  ==  t_stop 
break 

end 

if  t_stepvalue  >  t_stop 

t_step (count+1 )  =  t_stop; 
h (count  +1)  =  t_stop  -  t_step (count) ; 
h2  (count2  +1)  =  t_stop  -  t_step (count)  ; 
t_stepcount  =  t_step(count  +1); 
t_stepcount  =  t_step(count  +1); 
delt  =  h (count+1) ; 

[rl,cl]  =  size  (t  step) ; 
for  ii  =  l:cl 

t_diff(ii)  =  t_step(end)  -  t_step(ii) 

end 

xcoarse  =  x_coarse (1 : count) ; 
x  half  =  x  half (1 : count) ; 
x  fine  =  x  fine (1 : count) ; 

end 

end 

toe 

double_count 
half  count 

total_count  =  count2 

%  The  Exact  Solution: 

delt2  =  .00001; 

k2  =  100  +  delk; 

wnChk  =  sqrt(k2/m); 

fn  =  wnChk/ (2*pi) ; 

wd  =  wnChk  *  sqrt(l  -  zA2); 

omega  =  4; 

W  =  2  *  pi  *  omega; 
f Amp  =  1.0; 
f  =  fAmp*sin (W*time) ; 

%  f  =  ones (size (time) ) ; 

hchk  =  fModalIRFkernel (k2,m, z, time) ; 

xchk  =  conv (hchk, f ) ; 
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xchk  =  xchk ( 1 : length (hchk) )  *  delt2; 

%  Plotting  the  results: 

%  Change  default  axes  fonts. 
set(0,  ' Def aultAxesFontName ' ,  ' Arial') 

set  (0,  ' DefaultAxesFontSize '  ,  14) 

%  Change  default  text  fonts. 
set(0,  ' Def aultTextFontname ' ,  'Arial') 
set (0,  ' DefaultTextFontSize '  ,  14) 
figure  ( 1 ) 

set (gcf , ' WindowStyle ' , 'Docked') 

plot(t  step ( 1 : end) , x  approx, ' dk ', time, xchk) 

legend ( ' x* (t) ' , ' Exact ' ) 

xlabel('Time  (sec) ') 

ylabel ( ' x* (t)  '  ) 

title('x*(t)  vs.  Time') 

grid  on 

figure  (2 ) 

set (gcf, ' WindowStyle ; , ; Docked ' ) 

[rl,cl]  =  size(h2); 

tt  =  l:cl; 

plot (tt, h2 , ' kd ' ) 

xlabel ( ' Step  Number') 

ylabel (' Step  Size') 

title (' Step  Size  vs.  Step  Number') 

grid  on 

figure ( 3 ) 

set (gcf, ' WindowStyle ; , ; Docked ' ) 

[rl,cl]  =  size(h2); 

tt  =  l:cl; 

plot (tt, h2 , ' kd ' ) 

xlabel (' Step  Number') 

ylabel (' Step  Size') 

title (' Step  Size  vs.  Step  Number') 

grid  on 

axis ( [0, 250,0, .025]) 

%  Determining  the  error: 

[rowl,coll]  =  size (t_step) ; 
for  rr  =  1 : coll 

xact  find  =  find(time>=t  step(rr)); 
xact  min  =  min (xact  find) ; 
xact  sol (rr)  =  xchk (xact  min); 

end 

max^diff  =  max (abs (xact  sol  -  x  approx)); 

error  eval  =  max (abs (xact^sol  -  x  approx) ) /max (abs (xact  sol)) 
error  eval2  =  min (abs (xact_sol  -  x  approx) ) /max (abs (xact  sol)) 


S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  function  is  the  coarse  solution  as  derived  by  Gordis  and  Neta 
%Created  by  Keenan  Coleman  on  8/25/14 

S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

function  [x^coarse, x^coarsevalue]  = 

x  course  solution (h, x^coarse, base  result, kern^elem, kernl , count, delk) 

lhsCoef  f__coarse  =  1  +  0 . 5*delk*h  (end)  *kernl  (end)  ; 
for  ss  =  1: count 
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h_term(ss)  =  (h(ss+l)  +  h(ss)); 

end 

x_product  =  0 . 5*delk*h_term. *kern  elem.*x  coarse; 
x  sum  =  base  result  -  sum (x_product) ; 
xcoarse (count+1 )  =  x_sum/lhsCoef f_coarse; 
x  coarsevalue  =  x_sum/lhsCoeff_coarse; 

end 


oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  function  is  the  half  solution  as  derived  by  Gordis  and  Neta% 
%Created  by  Keenan  Coleman  on  8/25/14  % 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

function  [x  half,x  halfvalue]  = 

x  half  solution (h,x  coarse, x  half , kernhalf, kernl , count, base  result  half 
, delk) 

lhsCoef f_half  =  1  +  0 . 25*delk*h (end) *kernl (end)  ; 
x_product  half  =  0; 

if  count  >=  2; 

for  ss  =  l:count-l 

h^term(ss)  =  (h(ss+l)  +  h(ss)); 

end 

x_product  half  =  0.5*delk*h  term. *kernhalf (1 : count  - 
1 ). *x_coarse ( 1 : count  -1)  ; 
end 

x  sum  =  -sum (x_product  half)  -  0 . 5*delk* (h (count)  +  0.5*h (count  + 

1 )) *kernhalf (count) *x_coarse (count)  +  base  result  half; 
x  half (count+1 )  =  x  sum/lhsCoeff  half; 
x  halfvalue  =  x  sum/lhsCoeff  half; 

end 


ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  is  the  half  solution  as  derived  by  Gordis  and  Neta% 

%Created  by  Keenan  Coleman  on  8/25/14  % 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

function  [x_fine,  x  finevalue]  = 

x  fine  solution (h,x  coarse, x  fine,x  halfvalue, base  result, kern  elem, ker 
nl ,  kernhalf, count, delk) 

lhsCoef f_f ine  =  1  -  0 . 25*h (end) *kernl (end) ; 
x_product  =  0; 
if  count  >=  2; 

for  ss  =  l:count-l 

h_term(ss)  =  (h(ss+l)  +  h(ss)); 

end 

x_product  =  0.5*delk*h  term. *kernl ( 1 : count  -1 ). *x_coarse ( 1 : count 

-Dr- 

end 

x  sum  =  -sum (x_product)  -  0 . 25*delk* (h (count+1 )  + 

2*h (count) ) *kern  elem (count) *x_coarse (count)  - 
0 . 5*delk*h (end) *kernhalf (count) *x  halfvalue  +  base  result; 
x  fine (count+1 )  =  x  sum/lhsCoeff  fine; 
x  finevalue  =  x  sum/lhsCoeff  fine; 
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end 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  function  calculates  the  baseline  response. % 

%Created  by  Keenan  Coleman  on  8/25/14  % 

9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'S-9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooo 

function  [  x  base  ]  =  x  baseline (t  start,  t  stop,  min_step,  m,  wn,  wd, 
z) 

%  Initializing  the  variables: 
t  =  t_start :min_step : t_stop; 

%  The  kernel  function: 

kern  value  =  exp (-z*wn* (t) )  .*  sin (wd  *  (t) )  /  (m*wd) ; 

%  Periodic  forcing  function  parameters: 
omega  =  4; 

W  =  2  *  pi  *  omega; 
f Amp  =  1.0; 

%  The  forcing  function  (this  is  the  only  place  the  forcing  function  is 
needed) : 

[force^value]  =  fAmp*sin (W*t) ; 

%  force_value  =  ones (size (t) ) ; 

%  The  convolution  result: 

x  base  =  convfkern  value,  force^value) ; 

x  base  =  x  base ( 1 : length ( kern  value))  *  min^step; 

end 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  function  calls  the  given  base  resonse  for  the  given  time-step 
%Created  by  Keenan  Coleman  on  8/25/14 

S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

function  [base_result, base_result_half ]  =  x_base_call (time_start, 
min_del, time_stop,  t_stepl, t_step2, x_base) 

xx  =  time  start :min  del: time  stop; 

cl  =  find(xx  >=  t_step2); 

c2  =  find(xx  >=  (t_stepl  +  t_step2)/2); 

cl  =  min (cl ) ; 

c2  =  min (c2 ) ; 

base_result  =  x_base(cl); 

base_result_half  =  x_base(c2); 

end 


ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  function  evaluates  the  kernel  at  each  time-step  and  half% 
%time-step  needed  for  the  algorithms  derived  by  Gordis  and  % 
%Neta .  % 

%Created  by  Keenan  Coleman  on  8/25/14  % 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

function  [  kern  elem, kernl , kernhalf  ]  = 
kernel_eval (m, wn, wd, z, t  diff , t_step, count) 

[rowl , columnl ]  =  size  (t_diff) ; 
for  ii  =  1: (columnl-1) 

t_half(ii)  =  (t_step (count) +t_step (count+1 )) /2  -  t_step(ii); 
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end 

kern  =  @ (t  diff ) (1/ (m*wd) ) *exp (-z*wn*t  dif f ) . *sin (wd*t  diff); 

kernl  =  kern(t  diff); 

kernhalf  =  kern(t  half); 

kern  elem  =  kernl (1 : end-1) ; 

end 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  function  evaluates  the  IRF  for  values  of  t% 
%Created  by  Keenan  Coleman  on  8/25/14  % 

9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'S-9'9'9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9-9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooo 

function  ModelRF  =  fModalIRFkernel (k,m, z, t,  tau) 


wn  =  sqrt (k/m) ; 

wd  =  wn  *  sqrt(l  -  zA2); 

if  nargin  ==  4 

ModelRF  =  exp (-z*wn* (t) )  . * 


elseif  nargin 
ModelRF  = 


==  5 

exp (-z*wn* (t-tau) 


End 


sin  (wd  *  (t)  )  /  (m*wd)  ; 

.*  sin (wd  *  (t-tau))  / 


(m*wd) 
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APPENDIX  G.  MATLAB  CODE  FOR  THE  FINITE  ELEMENT 
MODEL  OF  THE  ALUMINUM  CANTILEVERED  BEAM 


9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9-9'9'9'9'9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  is  the  finite-element  beam  code  that  gives  the  exact  solution 
%used  to  compare  the  results  of  all  of  the  MDOF  models. 

%Created  by  Keenan  Coleman  on  8/25/2014 

9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


close  all; 
clear  all; 
clc; 

format  long; 

%  Problem  givens 
modulus  =  10e6; 
density  =  2.5879e-4  ; 
beamLength  =  120; 
beamHeight  =  12; 
beamWidth  =  4; 
numelem  =  5; 
dampingratio  =  .02; 

rootdamp  =  ( 1-dampingratio A2 ) A  (  .  5 ) ; 
area  =  beamHeight*beamWidth; 
moment  =  beamWidth*beamHeightA3/12; 

1  =  beamLength/numelem; 
alpha  =  density  *  area; 
numdof  =  numelem*2+2; 
precisefreq  = 

3.52/ (2*pi) *sqrt (modulus*moment/ (beamLengthA4*area*density)  )  ; 
t  =  0:. 0005:. 5; 

%  This  applies  the  necessary  applied  loads  and  moments. 
forceVec  =  zeros (numdof, 1) ; 
count  =  1; 
query  =  ' y ' ; 

while  query  ==  ' y '  &  count  <=  numdof 
nodeNum  =  11; 
while  nodeNum  >  numdof 

fprintf('The  input  exceeds  vector  dimensions\n ' ) ; 
nodeNum  =  input ('Enter  the  applicable  DOF:  '); 

end 

forceVec (nodeNum)  =  1000; 
count  =  count  +1; 
if  count  <=  numdof 
query  =  ' n ' ; 
end 

end 

Pknot  =  forceVec (nodeNum, 1 ) ; 

M  =alpha* [13/35*1,  11/210*1A2,  9/70*1,  -13/420*1A2;  11/210*1A2, 
1/105*1A3,  13/420*1A2, -1/140*1A3;  9/70*1,  13/420*1A2,  13/35*1,  - 
11/210*1A2;  -13/420*1 A2,  -1/140*1A3,  -11/210*1A2,  1/105*1A3]; 
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K  =  modulus*moment* [12/1A3,  6/lA2,  -12/1A3,  6/lA2;  6/lA2,  4/1,  -6/lA2 
2/1;  -12/1A3,  -6/1 A2 ,  12/1A3,  -6/lA2;  6/lA2,  2/1,  -6/lA2,  4/1]; 


%  This  builds  the  global 
response . 

mGlobal  =  zeros (numdof ) ; 
kGlobal  =  zeros (numdof) ; 
s  —  0 ; 

for  j  =  l:numdof/2-l 
for  ii=l:4 

for  k=l : 4 

kGlobal ( ii+s , 
mGlobal ( ii+s , 
end 

end 
s=s+2 ; 

end 


stiffness  and  mass  matrices  for  the  baseline 


k+s ) 
k+s ) 


K(ii,  k)  +  kGlobal  ( ii  +  s ,  k+s); 
M(ii,  k)  +  mGlobal ( ii+s ,  k+s); 


%  This  builds  the  global  stiffness  matrix  for  the  modified  response. 
mGlobal2  =  zeros (numdof) ; 
kGlobal2  =  zeros (numdof) ; 
s  =  0  ; 

for  j  =  l:numdof/2-l 
for  ii=l:4 

for  k=l : 4 

if  j  ==  3  %  Modified  element  number. 

kGlobal2 ( ii+s ,  k+s)  =  1.5*K(ii,  k)  +  kGlobal2 ( ii+s , 

k+s)  ; 

else 

kGlobal2 ( ii+s ,  k+s)  =  K(ii,  k)  +  kGlobal2 ( ii+s ,  k+s); 

end 

end 

end 
s=s+2 ; 

end 


%  This  is  to  impose  boundary  conditions  on  the  global  mass  and 
stiffness 

%  matrices  by  restraining  the  applicable  mode. 

restrain  =  ' y ' ; 
while  restrain  ==  ' y' 
rNode  =  1; 
dofl  =  2*rNode  -1; 
dof2  =  2*rNode; 
kGlobal (dofl : dof 2 , :)=[]; 
kGlobal ( : , dofl : dof 2 )  =  [ ]  ; 
kGlobal2 (dofl : dof 2 ,  :)  =  []; 
kGlobal2 ( : , dofl : dof 2 )  =  [ ]  ; 
mGlobal (dofl : dof 2 , :)=[]; 
mGlobal ( : , dof 1 : dof 2 )  =  [ ]  ; 
f orceVec (dofl : dof 2 )  =  [  ]  ; 
restrain  =  ' n ' ; 

end 
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%  This  is  the  calculation  for  the  tip  displacement  and  tip  rotation 
using  modal  displacement. 

[phi, lam]  =  eig ( kGlobal , mGlobal ) ; 

[phi2,lam2]  =  eig ( kGlobal2 , mGlobal ) ; 

[row, col]  =  size (phi); 

%Tip  Displacement  for  a  step  load. 

%  Baseline  response 
c  =  1; 

[rowl,coll]=  size(t); 
x  =  zeros (row, coll ) ; 
for  j  j  =  1 : row 
for  ii  =  t 

qinit(jj,c)  =  exp (-dampingratio*sqrt (lam ( j j , j j ) ) *ii) * ( (- 
dampingratio*phi (nodeNum- 

2,jj)*Pknot/( sqrt ( lam ( j  j , j  j ) ) A2*rootdamp) ) *sin ( sqrt ( lam ( j  j , j  j ) ) *rootdam 
p*ii)  .  .  . 

- (phi (nodeNum- 

2 , j  j ) *Pknot/sqrt ( lam ( j  j , jj) ) A2) *cos ( sqrt ( lam ( j  j , j  j ) ) *rootdamp*ii ) ) +phi ( 

nodeNum-2 , j  j ) *Pknot/sqrt (lam(jj,jj))A2; 

c  =  c+1; 

end 

c=l  ; 

end 

%  Modified  Response 
count  =  1; 

[row2,col2]=  size(t); 
x2  =  zeros (row, col2 ) ; 
for  j  j  =  1 : row 
for  ii  =  t 

qinit2 (j j  ,  count)  =  exp (-dampingratio*sqrt (lam2 ( j j  ,  j j ) ) *ii) * ( (- 
dampingratio*phi2 (nodeNum- 

2,jj)*Pknot/( sqrt ( lam2 ( j  j , j  j ) ) A2*rootdamp) ) *sin ( sqrt ( lam2 ( j  j , j  j ) ) *rootd 
amp  * i i )  .  .  . 

- (phi2 (nodeNum- 

2 , j  j ) *Pknot/sqrt ( lam2 (jj,jj))A2)*cos( sqrt ( lam2 ( j  j , j  j ) ) *rootdamp*ii ) ) +ph 
i2 (nodeNum-2 , j  j ) *Pknot/sqrt ( lam2 ( j  j  ,  j  j ) ) A2  ; 
count  =  count+1; 
end 

count=l ; 
end 


%%  Creating  the  vectors  to  be  plotted 
xplot  =  0; 
xplot2  =  0; 
for  ii  =  1 : row 

xplot  =  xplot  +  qinit (ii,  : ) *phi (3  ,ii); 
xplot2  =  xplot2  +  qinit2 ( ii , : ) *phi2 ( 3  , ii); 
end 

%  Plot  these  values  against  t  to  get  a  plot  of  the  system  response. 
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APPENDIX  H.  MATLAB  CODE  FOR  A  NON-ADAPTIVE,  NON¬ 
DERIVATIVE  APPROXIMATION  OF  THE  ALUMINUM 
CANTILEVERED  BEAM  (SUPPORTING  FUNCTIONS  INCLUDED) 


S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  is  the  coarse  solution  to  the  MDOF  beam  problem  that  exhibited% 
%conditional  stability  for  varying  stiffness  modifications  and  time  % 
%steps.  % 

%Created  by  Keenan  Coleman  8/25/14  % 

S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%  Initializing  the  time  variables: 

t_start  =  0; 

t_end  =  . 5 ; 

de 1 t  =  0.01; 

deltl  =  0.00001; 

min_step  =  0.00001; 

timel  =  t  start : deltl : t  end; 

t  min  =  t  start :min  step:t_end; 

nStep  =  length (t  min); 

nStep2  =  length (timel ) ; 

max  error  =  0.0001; 

%  Initializing  the  initial  conditions: 
x  initial  =  zeros  (4,1) ; 
h_step  =  [0  (delt  -  t_start) ] ; 
x  coarse  =  [x  initial]; 
x  half  =  [x  initial]; 
x  fine  =  [x  initial]; 
x  approx  =  [x_initial]; 
base_result_plot  =  [x_initial] ; 
count  =  1; 

t_step  =  [0  delt] ; 
tdiff  =  [delt  0]; 
t_stepcount  =  t_start  +  delt; 
count  =  1; 

%  Initializing  the  physical  parameters: 
wn  =  sqrt (diag ( lam) ) ; 
fn  =  wn/2/pi; 

wd  =  wn  . *  sqrt(l  -  dampingratio . A2 ) ; 
z  =  dampingratio*ones (row,  1 )  ; 
dk  =  0.5*K; 

%  Initializing  the  forcing  function: 
f (row,  :)=  zeros ( 1 , nStep) ; 
f ( (row  -  1),:)  =  1000*ones ( 1 , nStep) ; 

%  Calculate  IRF  matrix  for  baseline  system  (no  mods) : 
H  =  fHmatrix  rev3 (wn, wd, phi, z, t  min) ; 
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%  Calculate  baseline  response  using  straight  convolution  with  H  matrix 
xBase  =  MDOF  base_rev3 (H, f , min  step,nStep); 

%  Initiating  the  coarse  solution  algorithm: 
while  t_stepcount  <=  t  end 

%  Calling  the  necessary  values  of  the  baseline  solution: 

[base_result, base_result_half , base_result_plot]  = 

MDOF  base  call  rev3 (t_min, 

t_step (count) , t_step (count+1 ) , xBase, base_result_plot, count) ; 

%  Evaluating  the  kernel  matrix  over  the  adapted  time-step: 

[h  matrix]  =  h  kernel_rev3 (wn, wd, phi, z, t  diff ) ; 

%  Implementing  the  coarse  solution: 

[  x^coarse, x^coarsevalue  ]  =  MD0F_coarse_rev3 (h  matrix, 
base_result, x^coarse, dk, h_step, count) ; 

%  Implementing  a  non-adaptive  time-step... 
count  =  count  +  1 
h_step (count+1 )  =  delt; 

t_step (count+1 )  =  t_step (count)  +  delt; 
t_stepcount  =  t_step (count+1 )  ; 
t_stepvalue  =  t_step (count)  +  delt; 

[rl,cl]  =  size (t_step) ; 
for  ii  =  l:cl 

t_diff(ii)  =  t_step(end)  -  t_step(ii); 

end 

end 

%  Change  default  axes  fonts. 

set ( 0  ,  ' Def aultAxesFontName ' ,  ' Arial ' ) 

set  (0,  ' DefaultAxesFontSize '  ,  14) 

%  Change  default  text  fonts. 

set(0,  ' Def aultTextFontname '  ,  'Arial') 

set (0,  ' DefaultTextFontSize '  ,  14) 

figure 

plot ( t_step ( 1 : end-1 ) , x^coarse (1, :) , ' r- . ' , t , xplot ,  k--',t,xplot2, ' b ' ) 
legend ( ' x* (t)  '  ,  ' Unmodified  (exact)  ' ,  'Modified  (exact)  ' ) 
xlabel('Time  (sec) ') 

ylabel ( ' x* (t)  (modified  displacement)') 
title('x*(t)  vs.  Time') 
grid  on 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  function  evaluates  the  IRF  matrix  multiplied  by  phi  and  phi 
%transpose 

%Created  by  Keenan  Coleman  on  8/25/2014 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


function  [  H  ]  =  f Hmatrix^rev3 (wn, wd, phi, z, t  min) 
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[rl,cl]  =  size(phi); 

ModalHmatrix  =  zeros (rl , rl , length (t_min) ) ; 

H  =  zeros (rl , rl , length (t  min)); 
for  ii  =  l:rl 

ModalHmatrix ( ii ,  ii ,: )  =  exp (-z (ii) *wn (ii) *t  min)  .*  sin(wd(ii)  * 

t  min)  /  wd (ii) ; 

end 

for  iStep  =  1  :  length (t  min) 

H(:,:,iStep)  =  phi  *  ModalHmatrix (:,:, iStep)  *  phi'; 

end 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%  This  function  evaluates  the  kernel  over  the  time  difference  vector% 
%Created  by  Keenan  Coleman  on  8/25/2014  % 

S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

function  [h  matrix]  =  h  kernel_rev3 (wn, wd, phi, z, t_diff ) 
phil  =  phi (3:6, : ) ; 

[rl,c2]  =  size  (phil); 

ModalHmatrix  =  zeros (c2 , c2 , length (t_diff) )  ; 
h  matrix  =  zeros (rl , rl , length (t  diff ) ) ; 

%  The  kernel  evaluated  at  the  actual  time-steps: 
for  ii  =  1 : c2 

ModalHmatrix ( ii ,  ii ,: )  =  exp (-z (ii) *wn (ii) *t  diff)  .*  sin(wd(ii)  * 

t_diff)  /  wd(ii); 

end 

for  iStep  =  1  :  length (t  diff) 

h  matrix (:,:, iStep)  =  phil  *  ModalHmatrix (:,:, iStep)  *  phil'; 

end 

9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  function  calls  the  baseline  response  value  at  the  current% 
%time-step.  This  value  has  been  calculated  outside  of  the  loop.% 
%Created  by  Keenan  Coleman  on  8/25/2014.  % 

9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'S-9'9'9'S-9'9'9'9'9'9'9'S-9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


function  [base_result, base_result_half , base_result_plot]  = 
MDOF  base  call  rev3 (t  min, 

t_stepl , t_step2 , xBase, base_result_plot,  count) 
cl  =  find(t_min  >=  t_step2); 
c2  =  find(t_min  >=  (t_stepl  +  t_step2)/2); 
cl  =  min (cl ) ; 
c2  =  min (c2 ) ; 

base_result  =  xBase(3:6,cl); 
base_result_plot ( : , count+1 )  =  base_result; 
base_result_half  =  xBase (:, c2 ); 
end 


ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  function  calculates  the  baseline  response  using  the  convolution . % 
%Values  will  be  called  as  necessary  in  the  loop.  % 

9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


function  [  xBase  ]  =  MDOF  base_rev3 (H, f , min_step, nStep) 
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[r2,c2] 
xBase  = 
for  i  = 
for 


end 


=  size (f ) ; 
zeros (r2 , nStep) ; 

1  :  r2 
j  =  1  :  r2 

tmp  =  conv (squeeze (H (i 
xBase  (i,  :)  =  xBase  (i,  : 


end 

end 


j  ,  : )  )  ,  f  ( j  ,  : )  )  ; 

+  tmp(l:nStep) 


*  min  step; 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  is  the  MDOF  coarse  solution.  This  is  set  up  to  go  adaptive. 
%Created  by  Keenan  Coleman  on  8/25/2014. 

S-S-S-S-S-S-S-S-S-S-S-S-S'S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


function  [  x^coarse, x^coarsevalue  ]  =  MDOF  coarse_rev3 (h  matrix, 
base_result, x^coarse, dk, h_step,  count) 

[rl,cl]  =  size (h  matrix ( :  ,  : , 1 ) ) ; 
x  coarse^sum  =  zeros (rl,l); 
for  ii  =  1: count 

x  coarse  sum  =  x_coarse  sum  +  . 5* (h_step ( ii )  + 
h  step  (ii+1 )  )  *h  matrix  (  :  ,  :  ,  ii )  *dk*x__coarse  (  :  ,  ii )  ; 
end 

x_coarsevalue  =  base_result  -  x_coarse_sum; 

x_coarse ( : , count+1 )  =  x_coarsevalue; 

end 
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APPENDIX  I.  MATLAB  CODE  FOR  A  NON-ADAPTIVE, 
DERIVATIVE  APPROXIMATION  OF  THE  ALUMINUM 
CANTILEVERED  BEAM  (SUPPORTING  FUNCTIONS  INCLUDED) 


ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  is  the  non-adaptive  derivative  approximation  of  the  aluminum 
%cantilevered  beam. 

%Created  by  Keenan  Coleman  on  8/25/2014 

S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%  Initializing  the  time  variables: 

t_start  =  0; 

t_end  =  .  1 ; 

delt  =  0.00001; 

min_step  =  0.00001; 

t  min  =  t  start :min  steprt  end; 

nStep  =  length (t  min); 

max  error  =  0.0001; 

%  Initializing  the  initial  conditions: 
x  initial  =  zeros  (4,1); 
h_step  =  [0  (delt  -  t_start) ] ; 
xdot_coarse  =  [x_initial] ; 
x  vec  =  [x  initial]; 
xapprox^sum  =  [x  initial]; 
base_result_plot  =  [x_initial] ; 
trap_base  =  [x_initial] ; 
count  =  1; 

t_step  =  [0  delt] ; 
tdiff  =  [delt  0]; 
t_stepcount  =  t_start  +  delt; 
count  =  1; 

%  Initializing  the  physical  parameters: 
wn  =  sqrt (diag ( lam) ) ; 
fn  =  wn/2/pi; 

wd  =  wn  .*  sqrt(l  -  dampingratio . A2 ) ; 
z  =  dampingratio*ones (row,  1 )  ; 
de 1 k  =  K  ; 
dk  =  delk; 

%  Initializing  the  step-forcing  function: 
f(10,:)=  zeros ( 1 , nStep) ; 
f  ( 9 ,  : )  =  1000*ones (1, nStep) ; 
fl  =  zeros (row, 1) ; 
fl  (9, 1)  =  1000; 

tic 

t  diff2  =  fliplr(t  min)  ; 

[IRF]  =  MDOF  trap_f IRF (wn,  z,  phi,  t  diff2); 
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[hdot  matrix]  =  vec  to  matrix (IRF,t  diff2,  phi); 

[indexl,  index2]  =  size (t_dif f2 ) ; 

%  Initiating  the  adaptive  algoritms  from  Dr.  Neta's  paper: 
while  t_stepcount  <=  t  end 

%  Evaluating  the  exact  baseline  value  at  the  current  step: 

[trap_base, trap_base_value]  = 

base_exact (wn, wd, phi , z , f 1 , trap_base, count, t_stepcount ) ; 

%  Evaluating  the  kernel  matrix  over  the  adapted  time-step: 
hdot  matrixl  =  hdot  matrix (:,:, index2-count : end) ; 

%  Creating  the  x  vector  and  approximation  values  required  for  the 
%  coarse  solution 
if  count  >=2 

[x^vec,  xapprox^vec, xapprox_sum]  =  MDOF  x_vector (count,  x  initial, 
h_step,  x_vec,  xdot_coarse, xapprox_sum) ; 
end 

%  Implementing  the  coarse  solution: 

[  xdot_coarse, xdot  coarsevalue  ]  =  MDOF_coarse  dot (hdot  matrixl 
, trap  base  value, xdot_coarse, x_vec, dk, h^step, count,  x  initial , xapprox^su 
m)  ; 

%  Implementing  a  non-adaptive  time-step... 
count  =  count  +  1 
h_step (count+1 )  =  delt; 

t_round  =  round (10e5* (t_step (count)  +  delt) ) /10e5; 
t_step (count+1 )  =  t_round; 
t_stepcount  =  t_step (count+1 ) ; 
t_stepvalue  =  t_step (count)  +  delt; 


end 

toe 

%  Plotting  the  approximated  solution  to  the  exact  solution  and  the 
baseline 
%  solution: 
figure 

%  Change  default  axes  fonts. 

set ( 0 , ' Def aultAxesFontName ' ,  ' Arial ' ) 

set (0, ' DefaultAxesFontSize ' ,  14) 

%  Change  default  text  fonts. 

set(0,  ' Def aultTextFontname ' ,  'Arial') 

set  (0,  ' DefaultTextFontSize '  ,  14) 

set (gef , ' WindowStyle ' , 'Docked') 

plot(t,xplot2,  ' k ' , t_step ( 1 : end- 2 ) , x_vec (1 ,  : )  ,  ' r-- ' ) 
title (' Non-adaptive  MDOF  (Derivative  Approach) (x* (t) ) ') 
legend (' Modified  (Exact) ', 'Approximation') 
xlabel('Time  (sec) ') 
ylabel ( ' x*  (t)  '  ) 
grid  on 


S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  function  calculates  the  exact  baseline  response  at  the% 
%current  time-step.  % 

%Created  by  Keenan  Coleman  on  8/25/2014  % 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 
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function  [trap_base, trap_base_value]  = 

base_exact (wn, wd, phi , z , f 1 , trap_base, count, t_stepcount ) 

[rl,c2]  =  size  (phi); 

ModalHmatrixl  =  zeros (rl , rl ) ; 

%  The  baseline  evaluated  at  the  actual  time-steps: 
for  ii  =  l:rl 

ModalHmatrixl ( ii , ii )  =  exp (-z (ii) *wn (ii) *t  stepcount)  .*  sin(wd(ii)  * 

t_stepcount)  /  wd(ii); 

end 

hdot  matrixll  =  phi  *  ModalHmatrixl  *  phi'*fl; 

hdot  matrixll ( 1 : 2 )  =  []; 

hdotmatrixll (5 : 8)  =  []; 

trap  base  value  =  hdot  matrixll; 

trap_base ( : , count+1 )  =  trap_base_value; 


oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  function  calculates  the  IRF  for  a  given  mode  over  a  % 
%time  interval.  The  structure  of  this  function  was  % 

%provided  by  Dr.  Joshua  Gordis  and  modified  by  Keenan  Coleman% 

S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


function  [ModelRF]  =  fModalIRF_MDOF_trap  F (wn,  zeta,  t) 
wd  =  wn  *  sqrt(l  -  zetaA2); 

ModelRF  =  (wd*cos (t*wd)  -  wn*zeta*sin (t*wd) ) . / (exp (t*wn*zeta) *wd) 


%  Used  with  permission 

function  [IRF]  =  MDOF_trap_f IRF (Wn,  Zeta,  phi,  time); 

Q, 

O 

%  Usage:  [IRF]  =  fIRF(Wn,  Zeta,  phi,  time); 

o, 

o 

%  Creates  matrix  whose  ROWS  are  the  IRF  constructed  from 
%  modal  parameters  passed  into  the  function. 

Q, 

O 

%  Function  uses  all  modes  passed  in,  and  will  generate  all  IRF  for 
unique 

%  (symmetric)  input-output  pairs  for  all  rows  in  [phi] . 

%  IRF  are  stored  in  "symmetric  column  storage"  (See 
f SymmetricStore .m) 


Wn : 


%  Zeta: 
modes ) 

o, 

o 

%  phi : 
coordinates 


Vector  of  natural  frequencies  (rad/sec) 

If  Wn(i)  <  0.1,  rigid  body  mode  is  assumed. 

Scalar,  or  vector  damping  ratio  (scalar  applied  to  all 


Mass  normalized  modal  matrix.  Num  rows  =  number  of 

Num  cols  =  number  of  mode 


to  be  used. 

o, 

o 

%  time:  Time.  Row  vector  of  sampling  points  for  IRF  evaluation. 

O, 

o 

%  Written  by  Joshua  H.  Gordis,  Ph.D. 

%  Rev  #2  -  zeta  can  now  be  a  vector  1/18/06 


81 


ndof 

nmodes 

nsymcol 

columns 

IRF 

matrix 
mode IRF 
isymcols 
being  nsymcol 


=  size (phi, 1) ; 

=  size (phi, 2) ; 

=  ndof  *  (ndof  +  1)  /  2; 

=  zeros (nsymcol, length (time) ) ; 

=  zeros (1, length (time) ) ; 

=  0; 


if  length (Zeta)  ==  1; 

Zeta  =  Zeta  *  ones (nmodes,  1) ; 

end 


%  Number  of 
%  Initialize 

%  Will  end  up 


for  irows  =  1  :  ndof; 

for  icols  =  irows  :  ndof; 

isymcols  =  isymcols  +  1; 
for  imodes  =  1  :  nmodes; 

modelRF  =  fModallRF  MDOF  trap  F (Wn ( imodes )  ,  Zeta ( imodes ) , 

time) ; 


IRF ( isymcols ,: )  =  IRF ( isymcols ,: )  +... 

phi ( irows , imodes )  *  phi ( icols , imodes )  *  modelRF; 


end;  %  End 

icnt  modes 

end;  %  End 

icnt  col_dof 

end;  %  End 

icnt  row^dof 

%  End  function  fIRF.m 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  function  creates  and  builds  the  array  of  x  values  needed% 

%to  evaluate  the  coarse  solution  % 

%Created  by  Keenan  Coleman  on  8/25/2014  % 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

function  [x_vec,  xapprox^vec, xapprox  sum]  =  MDOF_x  vector (count, 
x  initial,  h  step,  x  vec,  xdot^coarse , xapprox  sum) 

[rl,cl]  =  size (xdot_coarse) ; 
xelem  sum  =  zeros (rl,l)  ; 

%  This  creates  the  vector  of  x  values  needed  in  the  coarse  solution. 
This 

%  is  also  the  solution  to  the  problem, 
if  count  ==2 

xelem_sum  =  0 . 5* (xdot_coarse (:, count-1 )  + 

xdot_coarse (:, count) ) *h_step (count)  +  x_initial; 
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else 


xelem_sum  =  x_vec (:, count  -1)  +  0 . 5* (xdot_coarse (:, count-1 )  + 
xdot_coarse ( : , count) ) *h_step (count)  ; 
end 

x  vec ( : , count)  =  xelem  sum; 

%  This  creates  a  vector  of  approximations  needed  in  the  coarse 
solution . 

xapprox^vec  =  0 . 5* (xdot_coarse (:, count-1 )  + 
xdot_coarse ( : , count) ) *h_step (count)  ; 
xapprox^sum  =  xapprox^sum  +  xapprox^vec; 


ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  function  evaluates  the  coarse  approximation  for  the  derivative% 
%approach.  % 

%Created  by  Keenan  Coleman  on  8/25/2014  % 

9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

function  [  xdot_coarse, xdot^coarsevalue  ]  = 

MDOF^coarse^dot (hdot_matrix, 

base  result, xdot^coarse, x_vec, dk, h  step, count, x  initial , xapprox^sum) 
[rl,cl]  =  size (hdot_matrix (:,:,!)); 
xdot^coarse^suml  =  zeros (rl,l); 

lhs_coeff  =  eye(rl)  +  0.25*h  step (end) A2*hdot  matrix (:,:, end) *dk; 
for  ii  =  1 : count 

xdot_coarse_suml  =  xdot_coarse_suml  +  . 5*(h  step(ii)  + 
h  step (ii+1 ) ) *hdot  matrix (:,:, ii ) *dk*x  vec(:,ii); 
end 

xdot_coarse_sum3  =  0.5*h  step (end) *hdot_matrix (:,:, end) *dk*xapprox_sum 
xdot^coarse  sum4  = 

0 . 25*h_step (end) A2*hdot  matrix ( : , : , end) *dk*xdot_coarse ( : , end)  - 
0.5*h  step (end) *hdot_matrix (;,:, end) *dk*x  initial; 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  function  takes  the  vector  of  IRF  and  creates  a% 
%lower  triangular  matrix.  % 

%Created  by  Keenan  Coleman  on  8/25/2014  % 

ooooooooooooooooooooooooooooooooooooooooooooooooooooo 

function  [hdot  matrix]  =  vec  to  matrix (IRF, t  diff,  phi) 
[rl,cl]  =  size (t_dif f ) ; 

[r2,c2]  =  size  (phi); 

hdot  matrix  =  zeros (r2 , c2 , cl ) ; 

for  ii  =  l:cl 

A  =  zeros (r2 ) ; 

ind  =  f ind (tril (ones  (10) ) )  ; 

A(ind)  =  IRF ( : , ii) ; 

A  =  A  +  tril (A, -1)  ’  ; 
hdot  matrix (:,:, ii )  =  A; 
end 

hdot  matrix (1:2,:,:)  =  [ ]  ; 
hdot_matrix (5:8,:,:)  =  [ ] ; 
hdotmatrix ( : , 1 : 2 , : )  =  [ ] ; 
hdot_matrix ( : , 5 : 8 ,  : )  =  [ ]  ; 
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APPENDIX  J.  MATLAB  CODE  FOR  THE  MATRIX  FORMULATION 
TO  THE  DERIVATIVE  OF  THE  INTEGRAL  FORMULATION  FOR 
TRANSIENT  STRUCTURAL  ANALYSIS  (SUPPORTING 
FUNCTIONS  INCLUDED) 


Many  of  the  functions  used  for  this  program  can  be  found  in  Appendix  I. 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  is  the  matrix  formulation  for  the  derivative  of  the  integral 
%formulation  for  transient  structural  synthesis. 

%Created  by  Keenan  Coleman  on  8/25/2014 

S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%  Initializing  the  time  variables: 

delt  =  0.00001; 

t_start  =  delt; 

t_end  =  . 1-delt; 

t_diff  =  t_start : delt : t_end; 

t  diff2  =  0:delt:t  end; 

nStep  =  length (t_diff) ; 

nStep2  =  length (t  diff2); 

%  Initializing  the  physical  parameters: 
wn  =  sqrt (diag ( lam) ) ; 
fn  =  wn/2/pi; 

wd  =  wn  .  *  sqrt(l  -  dampingratio . A2 ) ; 
z  =  dampingratio*ones (row,  1 )  ; 
de 1 k  =  K  ; 
dk  =  delk; 

%  Initializing  the  step-forcing  function: 
f(10,:)=  zeros ( 1 , nStep) ; 
f(9,:)  =  1000*ones (1 , nStep) ; 
fl  =  zeros (row, 1) ; 
fl  (9, 1)  =  1000; 

%  Calculating  the  array  of  matries  to  be  called  in  the  matrix 
formulation 

[IRF]  =  MDOF_trap_f IRF (wn,  z,  phi,  t  diff2); 

[hdot  matrix]  =  vec  to  matrix (IRF, t  diff2,  phi); 

[  hdot^dk  ]  =  hdotmultiplydk (  hdot_matrix, dk, nStep2  ); 
hdot_dk_zero  =  eye(4)  +  hdot_dk ( : , : , 1 ) * . 25*deltA2 ; 
hdot  dk  =  deltA2*hdot  dk; 
hdot_dk ( : ,  : , 1 )  =  0 . 5*hdot_dk ( : ,  :  ,  1 )  ; 

%  Building  the  matrix 
nCset=  4; 

A2  =  zeros (nCset*nStep2 , nCset*nStep2 ) ; 
elementnumber  =  [hdot  dk  zero] ; 
colcount  =  1:4; 
columnbuild  =  0; 
tic 
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for  intstep5  =  2:nStep2 

elementsum  =  zeros  (4,4); 
colcount  =  colcount  +  4; 
for  intstep6  =  l:intstep5 
if  intstep6  ==  intstep5 

elementsum  =  elementsum  +  0.5  *  hdot^dk ( : , :,intstep6) 
else 

elementsum  =  elementsum  +  hdot  dk ( : ,  :,intstep6); 

end 

end 

elementnumber (:, colcount)  =  elementsum; 
columnbuild  =  columnbuild  +  1 

end 
%  toe 

rows 2  =  -3:0; 
cols2  =  -3:0; 

RowCnt2  =  0; 

ColCnt2  =  0; 
matrixbuild  =  0; 

%  tic 

for  iStep2  =  1  :  nStep2; 

RowCnt2  =  RowCnt2  +  1; 
rows2  =  rows2  +  4; 
cols2  =  -3  :  0; 
cols3  =  rows2; 

for  jStep2  =  1  :  iStep2; 

ColCnt2  =  ColCnt2  +  1; 
cols2  =  cols2  +  4; 

[RowCnt2  ;ColCnt2]; 

A2 (rows2 , cols2 )  =  elementnumber (:, cols3 ) ; 
cols3  =  cols3  -  4; 


end 

matrixbuild  =  matrixbuild  +  1 

end 
%  toe 

%  Building  the  vector: 
cnt2  =  0; 
rowindex  =  [1:4] ; 
x  initial  =  zeros  (4,1); 
trap_base  =  [x_initial] ; 

%  tic 

for  intStep2  =  t_start : delt : . 1 

[trap_base, trap_base_value]  = 
base_exact (wn, wd, phi , z , f 1 , trap_base, count, intStep2 ) ; 
base_vec (rowindex,  1 )  =  trap_base_value; 
rowindex  =  rowindex  +  4; 

end 
%  toe 
bb  =  0  ; 

%  tic 

velocity  =  A2\base_vec; 
toe 
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%  Creating  and  plotting  the  velocity  vector: 

velocity2  =  [ 0 ] ; 

for  intstep3  =  lrnStep 

velocity2 ( 1 , intstep3+l )  =  velocity (4*intstep3-l, 1) ; 

end 

plot (t_dif f2 ,  velocity2) 

%  Plotting  the  displacement  vector: 


displacement  =  [0]  ; 
for  intstep4  =  2:nStep2 

displacement2  =  0.5*delt 
velocity2 (intstep4) ) ; 

displacement ( 1 , intstep4 ) 
displacement2 ; 
end 

%  Change  default  axes  fonts. 
set(0,  ' Def aultAxesFontName '  , 
set(0,  ' DefaultAxesFontSize '  , 
%  Change  default  text  fonts, 
set ( 0 ,  ' Def aultTextFontname ' , 
set(0,  ' Def aultTextFontSize  '  , 
plot  (t, xplot2 , ' -*k  , tdif f2 
title (' Matrix  Formulation  of 
Equation ' )  ; 

xlabel('Time  (seconds)'); 
ylabel ( ' x* (t)  '  ) 
legend ( ' Jerri ' , ' exact ' ) 
toe 


(velocity2 (intstep4-l)  + 

=  displacement (1, intstep4-l)  + 


' Arial ' ) 

14) 

'Arial ' ) 

14) 

displacement, '--r') 

the  Derivative  of  the  Synthesis 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This  function  multiplies  the  array  of  h  dot  matrices  by  the  elemental 
%stiffness  matrix 

%Created  by  Keenan  Coleman  on  8/25/2014 

9'9-9'9'9'9'9'9'9'9-9'9'9'9-9'9'9-9'9-9'9'9'9'9-9'9-9'9'9-9-9-9'9'9'9'9'9'9-9-9'9-9-9'9-9'9-9'9'9'9'9-9'9'9-9'9'9'9-9'9-9'9-9-9'9'9'9-9'9'9-9- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

function  [  hdot_dk  ]  =  hdotmultiplydk (  hdot_matrix, dk, nStep  ) 
for  ii  =  1 : nStep 

hdot  dk(:,:,ii)  =  hdot  matrix ii ) *dk; 

end 
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